In [1]:

    

# from https://github.com/pymc-devs/pymc/wiki/Surplus


    

In [2]:

    

"""
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)


    

In [18]:

    

model = MCMC((P_0, Imean, y, MSP, EMSP, B1990, Imean, k, r, iq, isigma2, itau2))


    

In [19]:

    

model.sample(10000, thin=20, burn=100, burn_till_tuned=True)


    

    




 [-------------------------------------205%-------------------------------------] 30626 of 14900 complete in 47.6 sec

In [20]:

    

plot(model)


    

    




Plotting MSP
Plotting EMSP
Plotting r
Plotting B1990
Plotting k
Plotting iq
Plotting isigma2
Plotting itau2






    













    













    













    













    













    













    













    

In [21]:

    

import spacepy.plot as spp


    

In [23]:

    

model.trace('B1990')[:]


    

    
Out[23]:





array([ 331.78037038,  331.95860817,  349.30886594,  345.74898493,
        358.35814941,  393.48316446,  425.29959133,  362.08721853,
        356.52508144,  371.68694282,  354.08970903,  345.29519792,
        304.06513783,  353.89471066,  357.40973544,  349.10277597,
        337.22253869,  312.80110748,  282.26373797,  293.8001022 ,
        296.85822661,  261.94830769,  290.96504536,  289.75176032,
        308.6060611 ,  334.37462687,  385.67044507,  390.63691404,
        346.87773311,  357.68620579,  349.82149115,  332.98034726,
        330.39278549,  387.3366436 ,  404.11170362,  359.05190109,
        315.37710645,  314.39322213,  280.78304899,  278.96253865,
        296.25531011,  311.20266381,  312.58529194,  313.63825677,
        317.64911287,  355.49974374,  324.07337307,  308.98174071,
        343.67953114,  377.03363197,  364.21417811,  403.72718439,
        363.59847768,  366.66695607,  356.11429874,  296.61041715,
        280.57176936,  298.77714164,  304.91680267,  275.73109999,
        304.51488475,  318.47391042,  315.19325892,  326.11176119,
        339.13762867,  305.58878495,  291.38287154,  316.44058403,
        334.9274748 ,  317.70450587,  340.76787421,  335.93698587,
        339.37158679,  322.43943634,  333.34380003,  302.82170025,
        278.34866489,  264.56843831,  302.68574392,  336.42256238,
        299.76835101,  323.47184043,  330.01021217,  294.52109187,
        338.75363698,  234.19720488,  256.78493855,  266.54446348,
        266.09394958,  291.10566412,  268.20390049,  262.37175311,
        274.89820978,  252.32874929,  249.22234979,  305.00748031,
        254.11552201,  260.94921587,  278.41513302,  251.73478994,
        218.26951446,  220.85692912,  231.00346247,  251.85030527,
        238.95113287,  239.25270064,  220.56028789,  229.17950957,
        213.98949178,  249.19574285,  299.65936931,  302.58525569,
        311.5121957 ,  302.6118109 ,  307.20198287,  295.19623637,
        267.28169801,  295.60185725,  297.29689991,  274.79334963,
        314.52335517,  319.09490237,  296.68589681,  263.06079535,
        258.71757276,  283.80288718,  296.68839217,  314.49503662,
        326.09561162,  304.8854214 ,  314.7425007 ,  317.56928487,
        326.76414609,  298.58012367,  329.0102676 ,  318.05539048,
        314.48306792,  286.73666871,  284.42698245,  264.15051614,
        264.07288042,  276.50008012,  264.86746964,  245.81532107,
        266.26178234,  246.01303279,  216.58605813,  225.15470785,
        238.48991448,  268.71443982,  313.91498394,  302.97104929,
        291.37338465,  326.55243196,  296.55008706,  303.17495638,
        314.71213171,  294.98201412,  293.60252937,  365.00031988,
        368.0288441 ,  378.09432777,  438.38092522,  487.83411908,
        482.04668009,  464.79336749,  369.86641824,  316.46399563,
        264.5247415 ,  263.14840114,  246.22843688,  231.30348895,
        241.16309111,  276.80377808,  299.99400644,  332.37612497,
        304.67097779,  337.49098106,  343.77845952,  341.67747339,
        348.24546521,  348.83785859,  334.12744333,  305.72225828,
        317.77008809,  334.19359728,  325.48593381,  261.81430316,
        280.636235  ,  313.77071921,  291.44406705,  334.27854944,
        354.75519802,  343.74665946,  364.90169142,  423.57072791,
        371.57191259,  376.91332379,  357.76478003,  396.65608986,
        376.3352718 ,  368.19421007,  418.59680012,  394.28917838,
        383.24741306,  354.95580088,  348.54363298,  353.79462869,
        315.77835587,  347.85158617,  383.2925489 ,  362.90710502,
        394.8459376 ,  358.46721575,  407.55254887,  391.69466573,
        336.27642274,  336.31937134,  355.86946303,  339.99488194,
        312.8878555 ,  292.0260548 ,  320.5494005 ,  228.69203389,
        232.49081618,  248.65788614,  248.63628254,  251.3556944 ,
        244.32278265,  271.62452235,  274.76323083,  272.11293985,
        275.02393107,  264.96480277,  305.61230216,  319.54795029,
        281.78722703,  243.76750792,  215.19422547,  288.49414347,
        239.33722948,  242.34578873,  224.51157418,  237.3062295 ,
        235.65513859,  280.21873913,  253.83785116,  271.01219248,
        324.42085036,  295.44514723,  270.10451563,  291.77264391,
        296.15004823,  237.38826442,  247.58752636,  286.35812907,
        264.2728194 ,  250.95750465,  277.19459902,  277.7558142 ,
        309.27858937,  296.63303481,  279.7791638 ,  303.21287976,
        288.94946507,  275.96362529,  291.71434357,  297.24811265,
        344.67118922,  347.05813361,  357.8733788 ,  339.01109888,
        311.57034343,  417.95132177,  372.71040333,  313.20579702,
        325.96901109,  326.95150888,  309.34008424,  355.34363623,
        396.58779238,  375.86242318,  400.78691118,  432.36076135,
        478.28900415,  493.31061825,  468.95951406,  464.24212038,
        423.27174541,  464.94671074,  355.1755315 ,  392.034579  ,
        393.95575715,  402.54160131,  374.94749902,  310.97048577,
        323.7948332 ,  234.2834908 ,  260.16488601,  244.16173534,
        267.44234747,  286.50701906,  302.30197931,  322.15901153,
        299.75111024,  258.56926764,  252.75319772,  275.10760732,
        296.80154603,  280.71550545,  264.32947181,  245.10460571,
        233.82637751,  262.98979618,  220.1308481 ,  238.68348593,
        265.44217604,  272.71321502,  269.79904783,  276.76387373,
        296.05121317,  341.84500105,  323.81579416,  297.23380964,
        330.88303397,  304.20838299,  295.69294021,  281.91858554,
        259.69449167,  308.82301578,  290.97977011,  289.04879749,
        341.47477031,  390.00721689,  416.0912182 ,  336.81026228,
        321.92652212,  336.92778739,  289.13183904,  281.67329109,
        298.6610649 ,  281.02044694,  277.19553565,  270.23943198,
        302.82270677,  294.06262131,  305.83984779,  330.28575154,
        338.16682332,  388.07445002,  363.27253307,  370.1811208 ,
        355.31624862,  386.19025724,  385.29005423,  355.02865792,
        359.83858208,  331.59406629,  348.40741199,  353.61480118,
        303.23133172,  292.07429196,  292.72462429,  308.58208778,
        311.8524818 ,  326.58335256,  363.20545425,  395.73699181,
        443.39261333,  403.49727384,  391.59916833,  397.68886079,
        480.55063921,  600.69794267,  578.39368503,  551.91482705,
        525.98960244,  517.5802831 ,  541.18203339,  544.56096669,
        518.88452718,  532.28684257,  534.66672366,  477.19053943,
        480.42626207,  481.13952125,  511.42017308,  552.00550635,
        482.56587491,  454.51534638,  516.59782295,  568.03249265,
        517.34526557,  484.95286906,  438.67634057,  499.84687252,
        462.25323201,  465.88848645,  384.48570208,  355.51751935,
        339.84169396,  358.66473962,  421.0062528 ,  430.48205937,
        478.6714116 ,  408.43851292,  366.24577516,  391.65783931,
        288.84069337,  277.29715404,  310.47447054,  285.92108797,
        288.4461119 ,  280.75734419,  290.40231176,  296.18372728,
        299.77802749,  323.05327003,  369.27192943,  325.35423727,
        314.52938971,  298.93306125,  332.28588102,  299.67662666,
        296.78945451,  285.35387672,  307.68184106,  326.61622231,
        304.64674229,  304.07969462,  283.22057312,  320.38986236,
        339.82291616,  324.37046959,  297.2230637 ,  310.98926403,
        289.74974308,  288.68671584,  287.24590762,  276.99898669,
        314.30057337,  297.06210783,  296.58397635,  317.37027049,
        328.40875309,  310.40287749,  309.57531139,  282.67023043,
        278.36164105,  275.96522988,  279.09134551,  283.70379405,
        244.60427055,  242.54156102,  242.27018417,  238.05210626,
        235.04865968,  224.52398075,  232.84768956,  253.77811544,
        237.55037361,  246.73540386,  222.17882457,  218.29580324,
        234.30146814,  231.57654692,  254.47786063,  206.04549913,
        252.13276092,  308.21174732,  382.64163774,  366.04602886,
        338.1469274 ,  312.71551754,  302.18461557,  321.32138592,
        349.45489891,  333.15546463,  329.37042338,  286.11821047,
        318.48147807,  325.87747659,  324.21025474,  355.80904383,
        342.457374  ,  404.16967221,  347.444383  ,  332.42223109,
        311.80154975,  320.65819615,  306.17382654,  339.41813665,
        377.46693048,  399.73181839,  381.90120797])

In [ ]: