LatentOccupancy

"""
latent_occupancy.py

Simple model demonstrating the estimation of occupancy, using latent variables. Suppose
a population of n sites, with some proportion pi being occupied. Each site is surveyed, 
yielding an array of counts, y:

y = [3, 0, 0, 2, 1, 0, 1, 0, ..., ]

This is a classic zero-inflated count problem, where more zeros appear in the data than would
be predicted by a simple Poisson model. We have, in fact, a mixture of models; one, conditional
on occupancy, with a poisson mean of theta, and another, conditional on absence, with mean zero. 
One way to tackle the problem is to model the latent state of 'occupancy' as a Bernoulli variable 
at each site, with some unknown probability:

z_i ~ Bern(pi)

These latent variables can then be used to generate an array of Poisson parameters:

t_i = theta (if z_i=1) or 0 (if z_i=0)

Hence, the likelihood is just:

y_i = Poisson(t_i)

(Note in this elementary model, we are ignoring the issue of imperfect detection.)

Created by Chris Fonnesbeck on 2008-07-28.
Copyright (c) 2008 University of Otago. All rights reserved.
"""

# Import statements
from numpy import random, array
from pymc import MCMC, Matplot, Beta, Bernoulli, Lambda, Poisson, Uniform, deterministic

# Sample size
n = 100
# True mean count, given occupancy
theta = 2
# True occupancy
pi = 0.4

# Simulate some data data
y = [(random.random()<pi) * random.poisson(theta) for i in range(n)]

# Estimated occupancy
p = Beta('p', 1, 1)

# Latent variable for occupancy
z = Bernoulli('z', p, value=array(y)>0, plot=False)

# Estimated mean count
theta_hat = Uniform('theta_hat', 0, 100, value=3)

@deterministic(plot=False)
def t(z=z, theta=theta_hat):
    """Per-site Poisson count parameter"""
    
    return z*theta

# Poisson likelihood
counts = Poisson('counts', t, value=y, observed=True)

def main():
    # Run the model
    M = MCMC([p, z, theta_hat, t, counts])
    M.sample(110000, 10000, verbose=2)
    
    Matplot.plot(M)

if __name__ == '__main__':
    main()