Random effects band recovery model in PyMC

recovery.py

from pymc import *
from numpy import array, resize, diag, exp, zeros, prod, concatenate
import pdb

# Data
Releases = array([12,12,10,14,17,21,9,30,25,41,45,79,97,80,62,80,74,74,46,45])
NRCperiods = len(Releases) 
recdata = resize((
8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,
0,7,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,
0,0,8,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,13,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
0,0,0,0,14,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,2,
0,0,0,0,0,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,10,
0,0,0,0,0,0,6,1,1,0,1,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,21,5,1,0,0,0,0,0,0,0,0,0,0,3,
0,0,0,0,0,0,0,0,22,0,1,0,0,0,0,0,0,0,0,0,2,
0,0,0,0,0,0,0,0,0,35,3,0,0,0,0,0,0,0,0,0,3,
0,0,0,0,0,0,0,0,0,0,43,0,0,0,0,0,0,0,0,0,2,
0,0,0,0,0,0,0,0,0,0,0,65,2,1,0,0,0,0,0,0,11,
0,0,0,0,0,0,0,0,0,0,0,0,66,2,1,0,2,0,0,0,26,
0,0,0,0,0,0,0,0,0,0,0,0,0,48,13,1,0,0,0,0,18,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,44,4,1,0,0,0,13,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,3,0,1,0,15,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,55,4,1,0,14,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,7,2,29,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,30,4,12,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,31,14
), (20,21))


# Priors on survival parameters
meanS = Normal('meanS', 0, 0.01, value=0.0)
lnsdS = Normal('lnsdS', 0, 0.01, value=0.0)
# Priors on capture parameters
meanp = Normal('meanp', 0, 0.01, value=1.0)
lnsdp = Normal('lnsdp', 0, 0.01, value=0.0)

@deterministic
def sdS(lnsdS=lnsdS):
    """sdS = exp(lnsdS)"""
    return exp(lnsdS)
    
@deterministic
def sdp(lnsdp=lnsdp):
    """sdp = exp(lnsdp)"""
    return exp(lnsdp)
    
# Hyper distributions survival and probability of recovery
Sdev = Normal('Sdev', 0, 1, value=zeros(NRCperiods))
pdev = Normal('pdev', 0, 1, value=zeros(NRCperiods))

@deterministic
def S(meanS=meanS, Sdev=Sdev, sdS=sdS):
    """
    Random effect model for survival.
    
    S = invlogit(meanS + Sdev*sdS)"""
    return invlogit(meanS + Sdev*sdS)

@deterministic
def p(meanp=meanp, pdev=pdev, sdp=sdp):
    """
    Random effect model for recovery.
    
    p = invlogit(meanp + pdev*sdp)"""
    return invlogit(meanp + pdev*sdp)

@deterministic
def q(p=p, S=S, plot=False):
    """Joint probabilities of capture and survival"""
    
    # Calculate the diagonal
    qmat = diag(p*S)
    
    # Probabilities above diagonal
    for i in range(NRCperiods):
        for j in range(i+1,NRCperiods):
            # Probability of capture in the current recapture period
            qmat[i,j] = p[j]*S[j]*prod([S[k] * (1.-p[k]) for k in range(i,j)])
            
    
    # Concatenate probabilities of an animal never being seen
    return array([concatenate((x,[1.-sum(x)])) for x in qmat])

@observed
def recoveries(value=recdata, q=q):
    """
    Data likelihood for recoveries in each period.
    
    recoveries ~ multinomial(N, q)"""
    N=Releases
    return sum([multinomial_like(value[i,i:], N[i], q[i,i:]) for i in range(NRCperiods)])