# R program for calculating the maximum likelihood estimates and CL for the number of sika deer in Sorachi subprefecture in Hokkaido.

# 空知振興局におけるエゾシカ個体数を最尤推定するための従来型 R プログラム：

# 経験ジェフェリーズ事前分布法や data cloning 法の推定値と比較するため。
#（今の場合は簡単なモデルなので，このような従来型の最適化手法で求めることも可能だが，複雑なモデルになると経験ジェフェリーズ事前分布法を用いるしかない。）

# Kohji Yamamura (2018) Appropriate use of Bayes (1763) estimation for obtaining maximum likelihood estimates. Japanese Journal of Conservation Ecology 23:39-56.

# 山村光司 (2018) ベイズ推定法の適切な活用について－エゾシカ個体数推定の例－. 保全生態学研究 23:39–56

# -------------------------------------------
# Function for calculating likelihood (Likelihood)

Likelihood <- function(X,I,K){

P <- length(I);N <- numeric(P);U <- numeric(P-1);N[1] <- abs(X[1]);sigma <- abs(X[2]);sigma2 <-sigma*sigma
var <- matrix(sigma2,nr=P-1,nc=P-1);var <- var+diag(sigma2,P-1)
for (t in 1:(P-1)) IL[t] <- log(I[t+1])
for (t in 1:(P-1)) N[t+1] <- max(1.21*N[t] - K[t], 1.0E-6)
for (t in 1:(P-1)) U[t] <- log(N[t+1]/N[1])
Likeli <- -dmvnorm(x=IL, mean=U, sigma=var, log=TRUE)
return(Likeli)
} # End of Likelihood function

# -------------------------------------------
# Function for performing estimation (Para_est)

Para_est <- function(I,K,XInit,method="BFGS",pvalue=0.05){
Z <- optim(XInit,Likelihood,NULL,I,K,hessian = TRUE, method=method)
N <- abs(Z$par[1]);sigma <- abs(Z$par[2]);SE <- sqrt(diag(solve(Z$hessian)))
qvalue <- -qnorm(pvalue/2);LLimit_N <- N-qvalue*SE[1];ULimit_N <- N+qvalue*SE[1];LLimit_sigma <- sigma-qvalue*SE[2];ULimit_sigma <- sigma+qvalue*SE[2]

result <- cbind(c(N=N,sigma=sigma),c(LLimit_N,LLimit_sigma),c(ULimit_N,ULimit_sigma))
colnames(result) <- c("Estimates",paste("Lower",deparse(100*pvalue/2),"%CL",sep=""),paste("Upper",deparse(100*pvalue/2),"%CL",sep=""))
title <- paste("< Maximum likelihood estimates and ",deparse(100*(1-pvalue)),"% confidence interval >",sep="")
cat("\n",title,"\n",sep="")
print (result)
} # End of Para_est function

# -------------------------------------------
## Performing calculation ##
library(mvtnorm)
K <- c(0.3576521224,0.3610508872,0.5986849977,0.6319396895,0.7294055376,0)
I <- c(1.0,0.926101,1.111491,1.034397,1.263537,1.057444)
XInit <- c(4, 0.05)
Para_est(I,K,XInit)

# -------------------------------------------
## Results of calculation ##
#
# < Maximum likelihood estimates and 95% confidence interval >
#        Estimates Lower2.5%CL Upper2.5%CL
# N     2.66549489  2.24840468   3.0825851
# sigma 0.08067536  0.03070366   0.1306471