#Observation error surplus production model for albacore tuna with lognormal errors
###CAUTION - This code does not compile!!!
#############################################################################

model ObsErrorTuna;
{
######distribution of I's#####
for (i in 1:N) { Imean[i] <- log(Q*P[i]);
                      I[i] ~ dlnorm(Imean[i],itau2);}

#Calculate (deterministic) biomass trajectory
###################
P[1] <- 1;
for (i in 2:N) {
    P[i] <- P[i-1] + r*P[i-1]*(1-P[i-1]) - k*C[i-1];
    }

#####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#######
# Add lognormal prior corresponding to 10% and 90% of k at 1/300 and 1/80
k ~ dlnorm(-5.042905,3.7603664)I(0.001,0.2)
K <- 1/k;

#######Prior on Q#######
iq ~ dgamma(0.001,0.001);
q <- 1/iq;
Q <- q*K;

#Prior on tau2
#################################
itau2   ~ dgamma(1.708603,0.008613854);
tau2   <- 1/itau2;

MSP <- r*K/4;
}

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)