Salamanders

#!/usr/bin/env python
# encoding: utf-8
"""
salamanders.py

Example of occupancy estimation with constant probabilities. Data is from
repeated observations of salamander habitat sites, and is taken from MacKenzie et al. 2006
(http://books.google.co.nz/books?id=RaCmF9PioCIC).

Created by Chris Fonnesbeck on 2008-05-29.
Distributed under the MIT license: http://www.opensource.org/licenses/mit-license.php
"""

from pymc import *
from numpy import array

# Occupancy data
salamanders = array([[0,0,0,1,1], [0,1,0,0,0], [0,1,0,0,0], [1,1,1,1,0], [0,0,1,0,0], [0,0,1,0,0], [0,0,1,0,0], [0,0,1,0,0], [0,0,1,0,0], [1,0,0,0,0], [0,0,1,1,1], [0,0,1,1,1], [1,0,0,1,1], [1,0,1,1,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,0,0], [0,0,0,1,0], [0,0,0,1,0], [0,0,0,0,1], [0,0,0,0,1]])

# Number of replicates
k = 5

# Prior on probability of detetion
p = Beta('p', alpha=1, beta=1, plot=False)

# Prior on probability of occupancy
psi = Beta('psi', alpha=1, beta=1)

# Probability of not being detected on k occasions
not_detected = lambda_deterministic('not_detected', lambda p=p: (1.0-p)**k)

# Probability of no detections, conditional on occupancy
pc = lambda_deterministic('pc', lambda psi=psi, nd=not_detected: psi*nd)

# Probability of observing a zero -- sum of conditional probabilities
p0 = lambda_deterministic('p0', lambda psi=psi, pc=pc: (1.0-psi) + pc)

@observed(dtype=int)
def occupancy(value=salamanders, p0=p0, psi=psi, p=p):
    """Likelihood of observed occupancy data"""
    
    # Initialise likelihood sum
    like = 0.0
    
    # Loop over observations
    for s in value:
        
        if not sum(s):
            # Likelihood of no detections
            
            like += bernoulli_like(1, p0)
        
        else:
            # Likelihood of one or more detections
            
            like += bernoulli_like(1, psi)
            like += binomial_like(sum(s), k, p)
    
    return like