%matplotlib inline
import matplotlib.pyplot as plt
#from data import *
import pandas as pd
import numpy as np
import pymc as pm
pm.__version__
'2.3.3'
Y = pd.read_csv('fake_microtus.csv', header=None)
mass = Y.pop(5)
Y.shape
(348, 5)
Nobs, K = Y.shape
Nz = 125
# Agumentation size
Nz = 125
# Augmented data
Y = np.r_[Y, np.zeros((Nz,K))]
# Observed individuals
sighted = Y.sum(axis=1)>0
sighted_masked = np.ma.masked_equal(sighted,False)
nobs = len(sighted)
# Unobserved individuals
no_sightings = np.array([i for i,x in enumerate(sighted) if not x])
# Individual masses with missing values
missing_mass = np.r_[np.array(mass), [0.0]*Nz]
#mass = np.ma.masked_equal([[x]*K for x in missing_mass], 0.0)
mass_masked = np.ma.masked_equal(missing_mass, 0.0)
from pymc import Uniform, Lambda, Normal, Bernoulli, invlogit, MCMC, Matplot, Binomial, rbernoulli, potential, bernoulli_like
def generate_model():
# Covariate priors
a0 = Uniform('a0', lower=-10, upper=10, value=0)
a1 = Uniform('a1', lower=-10, upper=10, value=0)
# Mass mean
mu = Uniform('mu', lower=-10, upper=10, value=0.3)
# Mass SD
sigma_0 = Uniform('sigma_0', lower=0, upper=100, value=1.2)
# Mass precision
tau_0 = Lambda('tau_0', lambda sd=sigma_0: sd**-2)
# Imputed masses for agumented group
iMass = Normal('iMass', mu=mu, tau=tau_0, value=mass_masked, observed=True)
# Detection model
phi = Lambda('phi', lambda a0=a0, a1=a1, mass=iMass: invlogit([a0 + a1*mass for k in range(K)]))
# P(presence) for augmented groups
psi = Uniform('psi', lower=0, upper=1, value=0.5)
# Occupancy state for agumented group
Z = Bernoulli('Z', psi, value=sighted_masked, observed=True)
# Z = Lambda('Z', lambda psi=psi: np.r_[np.ones(Nobs), rbernoulli(psi, size=Nz)])
# Z_obs = Bernoulli('Z_obs', psi, value=np.ones(Nobs), observed=True)
# Z_aug = [Bernoulli('Z_aug_%i' % i, psi, value=1) for i in range(Nz)]
# Z = Lambda('Z', lambda z=Z_aug: np.r_[np.ones(Nobs), z])
# Detection given presence
muY = Lambda('muY', lambda Z=Z, p=phi: np.transpose(Z*p))
# Likelihood
Yi = Bernoulli('Yi', p=muY, value=Y, observed=True)
# Posterior estimate for population size
N = Lambda('N', lambda Z=Z: Z.sum())
return(locals())
from pymc import AdaptiveMetropolis, BinaryMetropolis, Slicer
M = MCMC(generate_model())
[M.use_step_method(pm.StepMethods.TWalk, p) for p in M.psi, M.mu, M.sigma_0, M.a0, M.a1]
[None, None, None, None, None]
M.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 284.8 sec
M.N.summary()
N: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 371.99 7.098 0.677 [ 358. 385.] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 359.0 368.0 372.0 376.0 387.0
Matplot.plot(M.N)
Plotting N
Matplot.plot(M.psi)
Plotting psi_0