############################
#101218-19ベイズ本輪読会@北海道大学
#
#このコードは飯島が勝手に作成した物であり、間違っている可能性があります。
############################
#11.2のギブズサンプラー
library(R2WinBUGS)
F <- matrix(c(50, 22, 14, NA, 28, 21, NA, NA, 25), ncol=3, byrow=TRUE)
NF <- matrix(c(78, 78, 79, 79, NA, 117, 117, 117, NA, NA, 187, 188), ncol=4, byrow=TRUE)

Nfound <- c(314, 351, 375)
R <- rep(400, 3)
N.study <- 3


#初期値
inits <- function() {list(
                     S = rep(0.5, 3),
                     Lambda = rep(0.5, 3),
                     NF = NF
                    )
         }
#結果で見るパラメータの設定
para <- c("S", "Lambda", "NF")

data <- c("F", "Nfound", "R", "N.study")


#モデル式（recapture.bugというファイルにこのモデル式を保存する）
model {

##########
##step2
##########
##Definition of indices d and A
#The number of died birds (found and not found)
d[1] <- F[1, 1] + NF[1, 1]
d[2] <- sum(F[1:2, 2]) + sum(NF[1:2, 2])
d[3] <- sum(F[1:3, 3]) + sum(NF[1:3, 3])

#The number of alive birds at the start of year j
A[1] <- R[1]
A[2] <- A[1] - d[1] + R[2]
A[3] <- A[2] - d[2] + R[3]


for (s in 1:N.study) {
     S[s] ~ dbeta(alphaS[s], betaS[s])
     alphaS[s] <- alpha1 + A[s] - d[s]
     betaS[s] <- beta1 + d[s]
}
alpha1 <- 1 #noninformative
beta1 <- 1 #noninformative


##########
##step3
##########
for (l in 1:N.study) {
     Lambda[l] ~ dbeta(alphaL[l], betaL[l])
}
alphaL[1] <- alpha2 + sum(F[1, 1:3])
alphaL[2] <- alpha2 + sum(F[2, 2:3])
alphaL[3] <- alpha2 + F[3, 3]
betaL[1] <- beta2 + Nfound[1] - NF[1, 4]
betaL[2] <- beta2 + Nfound[2] - NF[2, 4]
betaL[3] <- beta2 + Nfound[3] - NF[3, 4]

alpha2 <- 1 #non informative
beta2 <- 1 #non informative

##########
##step4
##########
##Calculation of cell probability
#For not found birds
pNF1[1] <- (1 - S[1])*(1 - Lambda[1])/phiNF[1]
pNF1[2] <- S[1]*(1 - S[2])*(1 - Lambda[2])/phiNF[1]
pNF1[3] <- S[1]*S[2]*(1 - S[3])*(1 - Lambda[3])/phiNF[1]
pNF1[4] <- S[1]*S[2]*S[3]/phiNF[1]
phiNF[1] <- (1 - S[1])*(1 - Lambda[1]) + S[1]*(1 - S[2])*(1 - Lambda[2])
          + S[1]*S[2]*(1 - S[3])*(1 - Lambda[3]) + S[1]*S[2]*S[3]
pNF2[1] <- (1 - S[2])*(1 - Lambda[2])/phiNF[2]
pNF2[2] <- S[2]*(1 - S[3])*(1 - Lambda[3])/phiNF[2]
pNF2[3] <- S[2]*S[3]/phiNF[2]
phiNF[2] <- (1 - S[2])*(1 - Lambda[2]) + S[2]*(1 - S[3])*(1 - Lambda[3])
          + S[2]*S[3]
pNF3[1] <- (1 - S[3])*(1 - Lambda[3])/phiNF[3]
pNF3[2] <- S[3]/phiNF[3]
phiNF[3] <- (1 - S[3])*(1 - Lambda[3]) + S[3]

#For found birds
pF1[1] <- (1 - S[1])*Lambda[1]/phiF[1]
pF1[2] <- S[1]*(1 - S[2])*Lambda[2]/phiF[1]
pF1[3] <- S[1]*S[2]*(1 - S[3])*Lambda[3]/phiF[1]
phiF[1] <- (1 - S[1])*Lambda[1] + S[1]*(1 - S[2])*Lambda[2]
          + S[1]*S[2]*(1 - S[3])*Lambda[3]
pF2[1] <- (1 - S[2])*Lambda[2]/phiF[2]
pF2[2] <- S[2]*(1 - S[3])*Lambda[3]/phiF[2]
phiF[2] <- (1 - S[2])*Lambda[2] + S[2]*(1 - S[3])*Lambda[3]
pF3 <- (1 - S[3])*Lambda[3] #only one value
i3found <- pF3 * R[3]

#Estimation of the number of not found birds (died and survived)
NF[1, 1:4] ~ dmulti(pNF1[], Nfound[1])
NF[2, 2:4] ~ dmulti(pNF2[], Nfound[2])
NF[3, 3:4] ~ dmulti(pNF3[], Nfound[3])

##Fitting for the number of found birds (data)
#The definition of the number of found birds
Found[1] <- R[1] - Nfound[1]
Found[2] <- R[2] - Nfound[2]
#Found[3] <- R[3] - Nfound[3]

#Fitting for data
F[1, 1:3] ~ dmulti(pF1[], Found[1])
F[2, 2:3] ~ dmulti(pF2[], Found[2])
F[3, 3] ~ dnorm(i3found, 1.0E-3)

} #model end



#計算
res <- bugs(data, inits, para, "recapture.bug",
            n.iter=10000, n.burnin=5000, n.thin=10, n.chain=3,
            working.directory=getwd(),
	    bugs.directory="C:/WinBUGS14/",
            clearWD = TRUE,
            debug = TRUE
       )