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

library(R2WinBUGS) # This should be replaced by library(R2OpenBUGS) in using OpenBUGS.

# ----------------------------------
# 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

# Observed population index
for (t in 1:P-1){
IL[t] <- log(I[t+1])
}

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

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

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

# Observation model
for (t in 1:P-1) { 
U[t] <- log(N[t+1]/N[1])
 }
IL[1:P-1] ~ dmnorm(U[1:P-1], ivar[1:P-1,1:P-1])
}

",fill=TRUE)
sink()

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

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

# ----------------------------------
# Calculation by WinBUGS (or OpenBUGS)
out <- bugs(data=sika.data, inits=inits, parameters.to.save=param, model.file="model.txt", n.thin=1, n.chains=1, n.burnin=10000, n.iter=110000, debug=FALSE, working.directory=getwd())
print(out, digits=5)

# ----------------------------------
# Calculation of skewness and quantiles
library(moments)
Skewness.out <- c(N2007=skewness(out$sims.list$N1trans),sigma=skewness(out$sims.list$Strans))
U95.out <- c(quantile(out$sims.list$N[,1],c(.975)),quantile(out$sims.list$sigma,c(.975)))
Median.out <- c(median(out$sims.list$N[,1]),median(out$sims.list$sigma))
L95.out <- c(quantile(out$sims.list$N[,1],c(.025)),quantile(out$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$sims.list$N1trans, col="gray",main = "", xlab="Transformed N2007", ylab="Density", cex.lab=1.3, cex.axis=1.2)
abline(v=median(out$sims.list$N1trans), lwd=3, col="red")
hist(out$sims.list$Strans, col="gray",  main = "", xlab="Transformed sigma", ylab="Density", cex.lab=1.3, cex.axis=1.2)
abline(v=median(out$sims.list$Strans), lwd=3, col="red")