PyMC model of koala sighting

koala.py

# ====================================================================================
# = Koala sightings example from Link and Barker 2009. Original WinBUGS model below. =
# ====================================================================================

import numpy as np
from numpy.ma import masked_array
from pymc import *

# Prior model weights, derived from Bayes factors   
pi = (0.016656, 0.052853, 0.045682, 0.087367, 0.091975, 0.221741, 0.19819,
0.284835, 0.0003505, 0.0003505)

# Count data from 3 observers, flattened from a 2x2x2 matrix to a tuple
# The last element represents number not seen by any observer.
count_data = np.array([43,14,12,8,7,7,6,-999])
count = masked_array(count_data, mask=count_data==-999)

# X-matrix: defines linear predictors for each count
X=np.reshape((1,1,1,1,1,1,1,
1,1,1,0,1,0,0,
1,1,0,1,0,1,0,
1,1,0,0,0,0,0,
1,0,1,1,0,0,1,
1,0,1,0,0,0,0,
1,0,0,1,0,0,0,
1,0,0,0,0,0,0), (8,7))

# Defines which parameters of a particular model are constrained to be zero
Xmod=np.reshape((1,1,1,1,0,0,0,
1,1,1,1,1,0,0,
1,1,1,1,0,1,0,
1,1,1,1,0,0,1,
1,1,1,1,1,1,0,
1,1,1,1,1,0,1,
1,1,1,1,0,1,1,
1,1,1,1,1,1,1,
1,1,1,1,1,1,1,
1,1,1,1,0,0,0), (10,7))

# Defines which parameters of a particular model are constrained to be
# equal to other parameters
PIM=np.reshape((0,1,2,3,4,5,6,
0,1,2,3,4,5,6,
0,1,2,3,4,5,6,
0,1,2,3,4,5,6,
0,1,2,3,4,5,6,
0,1,2,3,4,5,6,
0,1,2,3,4,5,6,
0,1,2,3,4,5,6,
0,1,1,1,4,4,4,
0,1,2,2,2,5,6), (10,7))

# Current model
M = Categorical('M', p=pi)

# M indices
MI = Lambda('MI', lambda M=M: (np.arange(10)==M).astype(int))

# All log-linear parameters
beta = Normal('beta', mu=0, tau=0.1, value=np.ones(len(count)-1)*0.5)

# Index log-linear parameters
b = Lambda('b', lambda M=M, beta=beta: Xmod[M]*beta[PIM[M]])

# Mean counts
lam = Lambda('lam', lambda b=b: np.exp(np.dot(X,b)))

# Poisson likelihood with missing element
C = Poisson('C', mu=lam, value=count, observed=True)

# Total population size
N = Lambda('N', lambda C=C: np.sum(C))


"""
# Original WinBUGS model

model
{
	Model~dcat(pi[])
	for(i in 1:K)
	{
		MI[i]<-equals(Model,i)
	}
	for(i in 1:Ncells)
	{
		count[i]~dpois(lam[i])
		log(lam[i])<-X[i,1]*b[1]+X[i,2]*b[2]+X[i,3]*b[3]
		+X[i,4]*b[4]+X[i,5]*b[5]+X[i,6]*b[6]+X[i,7]*b[7]
	}
	for(i in 1:Ncells-1)
	{
		beta[i]~dnorm(0,0.1)
		b[i]<-Xmod[Model,i]*beta[PIM[Model,i]]
	}
	N<-sum(count[])
}


list(K=10, Ncells=8, pi=c(0.016656, 0.052853, 0.045682, 0.087367,
0.091975, 0.221741, 0.19819, 0.284835, 0.0003505,0.0003505),
count=c(43,14,12,8,7,7,6,NA),
X=structure(.Data=c(
1,1,1,1,1,1,1,
1,1,1,0,1,0,0,
1,1,0,1,0,1,0,
1,1,0,0,0,0,0,
1,0,1,1,0,0,1,
1,0,1,0,0,0,0,
1,0,0,1,0,0,0,
1,0,0,0,0,0,0),.Dim=c(8,7)),
Xmod=structure(.Data=c(1,1,1,1,0,0,0,
1,1,1,1,1,0,0,
1,1,1,1,0,1,0,
1,1,1,1,0,0,1,
1,1,1,1,1,1,0,
1,1,1,1,1,0,1,
1,1,1,1,0,1,1,
1,1,1,1,1,1,1,
1,1,1,1,1,1,1,
1,1,1,1,0,0,0),.Dim=c(10,7)),
PIM=structure(.Data=c(1,2,3,4,5,6,7,
1,2,3,4,5,6,7,
1,2,3,4,5,6,7,
1,2,3,4,5,6,7,
1,2,3,4,5,6,7,
1,2,3,4,5,6,7,
1,2,3,4,5,6,7,
1,2,3,4,5,6,7,
1,2,2,2,5,5,5,
1,2,2,2,5,6,7),.Dim=c(10,7)))
#K = no. models, Ncells = no. cells in the table

"""