# Import modules
from pymc import *
import numpy as np
from get_data import *
# Constants
JUVENILE, ADULT = 0, 1
STARTYEAR, ENDYEAR = 2001, 2009
YEARS = ENDYEAR - STARTYEAR + 1
# Recovery rates from incidental take paper
recovery = dict(zip(regions, (0.8625, 0.8242, 0.4069, 0.5837)))
# Strong red tide years
tide_years = (2003, 2005, 2006)
# Specify Current model
AGE = 'Adults'
# Containers for variables
U = []
p = []
pi = []
X = []
X_like = []
Z = []
undetermined = []
theta = []
# Red tide effect factor
a = Uniform('a', lower=0, upper=2)
# Loop over regions
for area in regions:
# Upper St. John never gets red tide
k = len(causes) - 1*(area == 'USJ')
# Undetermined deaths by latent cause modeled as multinomial
# with Dirichlet prior on probabilities
p_i = Dirichlet('p_%s' % area, theta=np.ones(k), value=np.ones(k-1)/k)
p.append(p_i)
# Proportions of mortality for region i
pi_i = Dirichlet('pi_%s' % area, theta=np.ones(k))
pi.append(pi_i)
# Loop over years
for year in range(STARTYEAR, ENDYEAR+1):
# Subset data by year and age
cdata = data[(data["area"]==area) & (data['year']==year) & (data['age']==AGE)]
# Undetermined mortality
undetermined_i = cdata['undetermined']
undetermined.append(undetermined_i)
# Upper St. John never gets red tide
if area == 'USJ':
# Remove red tide as a cause
Z_i = cdata[causes[0:4] + [causes[-1]]].values[0]
else:
Z_i = cdata[causes].values[0]
# Undetermined deaths by latent cause modeled as multinomial
# with Dirichlet prior on probabilities
U_i = Multinomial('U_%s%i' % (area,year), n=undetermined_i, p=p_i)
U.append(U_i)
# Total deaths by cause, corrected for recovery rate
X_i = Lambda('X_%s%i' % (area,year), lambda u=U_i, z=Z_i: ((z + u)/recovery[area]).astype(int))
X.append(X_i)
# A red tide year ...
if (year in tide_years) and (area == 'SW'):
@deterministic
def theta_i(p=pi_i, a=a):
"""Adjust proportions in red tide years"""
theta = p.copy()
theta[4] += a
return theta/(1+a)
theta_i.__name__ += '_%s%i' % (area,year)
theta.append(theta_i)
# Normal year ...
else:
# No adjustment to proportions
theta_i = pi_i
@potential
def X_like_i(x=X_i, p=theta_i):
"The likelihood of total deaths by cause"
return multinomial_like(x, n=np.sum(x), p=extend_dirichlet(p))
X_like.append(X_like_i)
try:
del adult_model
except:
pass
adult_model = MCMC(locals())
adult_model.sample(20000, burn=10000)
[****************100%******************] 20000 of 20000 complete
adult_model.a.summary()
a: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.772 0.105 0.005 [ 0.587 0.991] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.589 0.699 0.765 0.836 0.996
Matplot.plot(adult_model.a)
Plotting a
[p.summary() for p in adult_model.pi]
pi_ATL: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.582 0.017 0.002 [ 0.546 0.61 ] 0.009 0.003 0.0 [ 0.003 0.015] 0.01 0.004 0.0 [ 0.004 0.018] 0.078 0.01 0.001 [ 0.057 0.099] 0.11 0.004 0.0 [ 0.103 0.116] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.547 0.572 0.583 0.594 0.612 0.004 0.007 0.009 0.011 0.017 0.005 0.008 0.009 0.012 0.018 0.059 0.071 0.078 0.085 0.102 0.102 0.107 0.11 0.113 0.116 pi_USJ: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.596 0.12 0.011 [ 0.379 0.806] 0.05 0.035 0.002 [ 0.002 0.121] 0.043 0.027 0.002 [ 0.001 0.095] 0.214 0.093 0.009 [ 0.038 0.398] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.364 0.505 0.606 0.695 0.796 0.006 0.024 0.041 0.068 0.14 0.005 0.023 0.039 0.059 0.113 0.067 0.15 0.2 0.264 0.442 pi_NW: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.61 0.054 0.005 [ 0.501 0.705] 0.011 0.01 0.001 [ 0. 0.031] 0.033 0.02 0.001 [ 0.003 0.07 ] 0.157 0.05 0.004 [ 0.065 0.249] 0.032 0.018 0.001 [ 0.004 0.068] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.508 0.572 0.612 0.648 0.712 0.0 0.003 0.008 0.015 0.037 0.006 0.02 0.029 0.044 0.084 0.077 0.12 0.153 0.189 0.265 0.007 0.018 0.029 0.043 0.075 pi_SW: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.445 0.014 0.001 [ 0.42 0.472] 0.032 0.006 0.001 [ 0.025 0.046] 0.008 0.004 0.0 [ 0.002 0.015] 0.086 0.011 0.001 [ 0.065 0.108] 0.332 0.014 0.001 [ 0.305 0.36 ] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.418 0.436 0.445 0.453 0.471 0.024 0.027 0.03 0.033 0.046 0.003 0.006 0.008 0.01 0.017 0.065 0.079 0.086 0.094 0.108 0.302 0.322 0.332 0.343 0.358
[None, None, None, None]