# Example of R-code for the estimation of sika deer population in Sorachi subprefecture.
# Empirical Jeffreys prior is used to obtain ML estimate.
# jags is used via R2jags package.

library(R2jags)

# ----------------------------------
# Initial parameters

# Exponents of Box-Cox transformation (for initial parameters)
lambdaN <- -2
lambdaS <- -0.5

initN1 <- 20  # Enter an appropriate initial population
initsigma <- 0.01  # Enter an appropriate initial sigma

initN1trans <- ((initN1^lambdaN)-1)/lambdaN
initStrans <- ((initsigma^lambdaS)-1)/lambdaS
inits <- function(){
list(N1trans=initN1trans,Strans=initStrans)}

# ----------------------------------
# Create a text file for describing the model
sink("model.txt")
cat("

model{

# Exponents of Box-Cox transformation
# Change the exponents to yield zero skewness.
lambdaN <- -2
lambdaS <- -0.5

# Range of prior uniform distribution
minN <- 1
maxN <- 500
minS <- 0.01
maxS <- 100

# Box-Cox transformation
LowerN <- (pow(minN,lambdaN)-1)/lambdaN
UpperN <- (pow(maxN,lambdaN)-1)/lambdaN
N1trans ~ dunif(LowerN,UpperN)
N1 <- pow((lambdaN*N1trans+1),(1/lambdaN))

LowerS <- (pow(minS,lambdaS)-1)/lambdaS
UpperS <- (pow(maxS,lambdaS)-1)/lambdaS
Strans ~ dunif(LowerS,UpperS)
sigma <- pow((lambdaS*Strans+1),(1/lambdaS))
sigma2 <-sigma*sigma

# Variance-covariance matrix for observation errors
for (t1 in 1:length(IL)){
for (t2 in 1:length(IL)){
variance[t1, t2] <- sigma2 + equals(t1,t2)*sigma2
}
}
ivar[1:length(IL),1:length(IL)] <- inverse(variance[,])

# Initialization of time series
N[1] <- N1

# State model
for (t in 1:length(IL)) { 
N[t+1] <- max(1.21*N[t] - K[t], 1.0E-6)
 }

# Observation model
for (t in 1:length(IL)) { 
U[t] <- log(N[t+1]/N[1])
 }
IL[1:length(IL)] ~ dmnorm(U[1:length(IL)], ivar[1:length(IL),1:length(IL)]) 
}
",fill=TRUE)
sink()

# ----------------------------------
# Data
sika.data <- list(K=c(3.576521224,3.610508872,5.986849977,6.319396895,7.294055376,0),IL=log(c(0.926101,1.111491,1.034397,1.263537,1.057444)))

# ----------------------------------
# Parameters to be monitored
param <- c("N1trans","N","Strans","sigma")

# ----------------------------------
# Calculation by jags
out <- jags(sika.data, inits=inits, parameters.to.save=param, model.file="model.txt", n.thin=1, n.chains=1, n.burnin=10000, n.iter=110000, working.directory=getwd())
print(out)

# ----------------------------------
# Calculation of skewness and quantiles
library(moments)
Skewness.out <- c(N2007=skewness(out$BUGSoutput$sims.list$N1trans),sigma=skewness(out$BUGSoutput$sims.list$Strans))
U95.out <- c(quantile(out$BUGSoutput$sims.list$N[,1],c(.975)),quantile(out$BUGSoutput$sims.list$sigma,c(.975)))
Median.out <- c(median(out$BUGSoutput$sims.list$N[,1]),median(out$BUGSoutput$sims.list$sigma))
L95.out <- c(quantile(out$BUGSoutput$sims.list$N[,1],c(.025)),quantile(out$BUGSoutput$sims.list$sigma,c(.025)))
cbind(Skewness=Skewness.out, Estimate=Median.out, "2.5 %"=L95.out, "97.5 %"=U95.out)

# ----------------------------------
# Posterior distribution in transformed scale
hist(out$BUGSoutput$sims.list$N1trans, col="gray",main = "", xlab="Transformed N2007", ylab="Density", cex.lab=1.3, cex.axis=1.2)
abline(v=median(out$BUGSoutput$sims.list$N1trans), lwd=3, col="red")
hist(out$BUGSoutput$sims.list$Strans, col="gray",main = "", xlab="Transformed sigma",ylab="Density", cex.lab=1.3, cex.axis=1.2)
abline(v=median(out$BUGSoutput$sims.list$Strans), lwd=3, col="red")