Surplus

"""
Example from Meyer & Miller (1999) Can. J. Fish. Sci. 56:1078
Created by Aaron MacNeil on 29.07.09 
"""

""" WinBUGS model from Russel Millar's website:
(http://www.stat.auckland.ac.nz/~millar/Bayesian/Surtuna.bugs)

#State-space surplus production model for albacore tuna with lognormal errors
#############################################################################
model surplustuna {
######distribution of I's#####
for (i in 1:N) { Imean[i] <- log(Q*P[i]);
                 I[i] ~ dlnorm(Imean[i],itau2);
                 B[i] <- K*P[i] }
######distribution of P's######
Pmean[1] <- 0;
P[1] ~ dlnorm(Pmean[1],isigma2)I(0.5,2.0)
for (i in 2:N) {
    Pmean[i] <- log(max(P[i-1] + r*P[i-1]*(1-P[i-1]) - k*C[i-1],0.01));
    P[i] ~ dlnorm(Pmean[i],isigma2);
    }
#####Prior on r######
#lognormal prior corresponding to 10% and 90% of r at 0.13 and 0.48 
r~dlnorm(-1.38,3.845);
#####Prior on k######
#Lognormal prior corresponding to 10% and 90% of k at 1/300 and 1/80
k ~ dlnorm(-5.042905,3.7603664);
K <- 1/k;
MSP<-r*K/4;
#####Prior on Q#####
iq ~ dgamma(0.001,0.001);
q <- 1/iq;
Q <- q*K;
#######Priors on isigma2 and itau2#####
a0<-3.785518;  b0<-0.010223; 
c0<-1.708603; d0<-0.008613854;
isigma2 ~ dgamma(a0,b0);
itau2   ~ dgamma(c0,d0);
Sigma2<-1/isigma2; Tau2<-1/itau2;
#######Project biomass into future#####
nyrs<-10
for (i in N:(N+nyrs-1)) {C[i]<-19}
for (i in (N+1):(N+nyrs)) {
    Pmean[i] <- log(max(P[i-1] + r*P[i-1]*(1-P[i-1]) - k*C[i-1],0.01));
    P[i] ~ dlnorm(Pmean[i],isigma2);
    B[i] <- K*P[i] }

}

list(N=23, 
     C=c(15.9, 25.7, 28.5, 23.7, 25.0, 33.3, 28.2, 19.7, 17.5, 19.3, 21.6, 23.1,
         22.5, 22.5, 23.6, 29.1, 14.4, 13.2, 28.4, 34.6, 37.5, 25.9, 25.3),
     I=c(61.89, 78.98, 55.59, 44.61, 56.89, 38.27, 33.84, 36.13, 41.95, 36.63, 
         36.33, 38.82, 34.32, 37.64, 34.01, 32.16, 26.88, 36.61, 30.07, 30.75, 
         23.36, 22.36, 21.91))

list(P=c(0.99,0.98,0.96,0.94,0.92,0.90,0.88,0.86,0.84,0.82,0.80,0.78,0.76,0.74,
         0.72,0.70,0.68,0.66,0.64,0.62,0.60,0.58,0.56),r=0.8, k=0.005, iq=5, 
     isigma2=100, itau2=100)


"""


#------------------------------------------------------------------ MODULES
from pymc import invlogit, Lognormal, deterministic, stochastic, data, Uniform, MCMC, lognormal_like, Gamma, normal_like, poisson_like, potential, uniform_like, observed, Lambda, Normal, AdaptiveMetropolis
from pymc.Matplot import plot
from numpy import zeros, ones, array, mean, std, exp, sqrt, pi, log
import pdb
    
#------------------------------------------------------------------ DATA

catch = array([15.9, 25.7, 28.5, 23.7, 25.0, 33.3, 28.2, 19.7, 17.5, 19.3, 21.6, 23.1, 22.5, 22.5, 23.6, 29.1, 14.4, 13.2, 28.4, 34.6, 37.5, 25.9, 25.3])
cpue = array([61.89, 78.98, 55.59, 44.61, 56.89, 38.27, 33.84, 36.13, 41.95, 36.63, 36.33, 38.82, 34.32, 37.64, 34.01, 32.16, 26.88, 36.61, 30.07, 30.75, 23.36, 22.36, 21.91])
year = range(1967,1990)
P_inits = array([0.99,0.98,0.96,0.94,0.92,0.90,0.88,0.86,0.84,0.82,0.80,0.78,0.76,0.74,0.72,0.70,0.68,0.66,0.64,0.62,0.60,0.58,0.56])
nyears = len(year)

#------------------------------------------------------------------ PRIORS

# Prior distribution of K
k = Lognormal('k', -5.042905, 3.7603664, value=0.005)

# Prior distribution of r
r = Lognormal('r', -1.38, 3.845, value=0.8)

# Prior distribution for q
iq = Gamma('iq', 0.001, 0.001, value=5.)

# Hyperparameters
isigma2 = Gamma('isigma2', 3.785518, 0.010223, value=100.)
itau2 = Gamma('itau2', 1.708603, 0.008613854, value=100.)

#------------------------------------------------------------------ MODEL

# Lognormal distribution of P's
Pmean0 = 0.
P_0 = Lognormal('P_0', mu=Pmean0, tau=isigma2, trace=False, value=P_inits[0])
P = [P_0]

# Recursive step
for i in range(1,nyears):
    Pmean = Lambda("Pmean", lambda P=P[i-1], k=k, r=r: log(max(P+r*P*(1-P)-k*catch[i-1],0.01)))
    Pi = Lognormal('P_%i'%i, mu=Pmean, tau=isigma2, value=P_inits[i], trace=False)
    P.append(Pi)

# Distribution of I's
@deterministic
def Imean(q=1/iq, K=1/k, P=P):
    return array([log((q)*K*p) for p in P])

@observed
@stochastic(plot=False)
def y(value=cpue, mu=Imean, tau=itau2):
    return lognormal_like(value,mu,tau)

#-------------------- Other tidbits

# Parameters on regular scale
K=Lambda('K', lambda k=k: 1/k)
q=Lambda('q', lambda iq=iq: 1/iq)

# Variances on regular scale
Sigma2=Lambda('Sigma2', lambda isigma2=isigma2: 1/isigma2)
Tau2=Lambda('Tau2', lambda itau2=itau2: 1/itau2)

# Quantities of interest
MSP=Lambda('MSP', lambda r=r, K=K: r*K/4)
EMSP=Lambda('EMSP', lambda r=r, q=q: r/(2*q))
@deterministic
def P1990(P=P[nyears-1], r=r, C=catch[nyears-1], k=k):
    return P+r*P*(1-P)-k*C
B1990=Lambda('B1990', lambda P1990=P1990, K=K: P1990*K)