# Example of R-code for the estimation of sika deer population in Sorachi subprefecture.
# dclone is used to obtain ML estimate.

library(dclone)
set.seed(1234)

# Object describing the model

Sorachi.model <- function(){

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

N1 ~ dunif(minN, maxN)
sigma ~ dunif(minS, maxS)
sigma2 <-sigma*sigma

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

for (j in 1:k) { 

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

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

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

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

}
}
# End of model

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

# Parameters to be monitored
param <- c("N1","sigma")

# Initial parameters
inits <- function(){
list(N1=30, sigma=0.1)}

# dclone data
data <- list(K=dcdim(K), I=dcdim(I), P=6, k=1)
dcdata <- dclone(data, n.clones = 50, multiply = "k", unchanged = "P")

# For OpenBUGS
Sorachi.ob50 <- bugs.fit(dcdata, param, Sorachi.model, inits=inits, n.chains=3, n.thin=1, n.burnin=10000, n.iter=110000, DIC = FALSE, debug=FALSE, program = "openbugs", working.directory=getwd())
summary(Sorachi.ob50)
cbind(Estimate=coef(Sorachi.ob50), confint(Sorachi.ob50))