とりあえず,
いろいろ準備して,
時系列データ解析の BUGS コードを書いてみた.
これまウミガメのやつよりかなりマシなものになっているはず
……
model
{
for (t in 1:T.max) {
Y[t] ~ dpois(m[t])
log(m[t]) <- log.m[t]
log.m[t] ~ dnorm(lm[t], tau[1])
lm[t] <- (log(Effort[t]) + log.cpue[t] + season.effect[Season[t]] - m.season)
}
Tau.noninformative <- 1.0E-04
log.cpue[1] ~ dnorm(lcm[1], Tau.noninformative)
log.cpue[2] ~ dnorm(lcm[2], tau[2])
log.cpue[3] ~ dnorm(lcm[3], tau[2])
for (t in 4:T.max) {
log.cpue[t] ~ dnorm(lcm[t], tau[2])
}
lcm[1] <- 0.0
lcm[2] <- log.cpue[1]
lcm[3] <- 2 * log.cpue[2] - log.cpue[1]
for (t in 4:T.max) {
lcm[t] <- 2.5 * log.cpue[t - 1] - 2 * log.cpue[t - 2] + 0.5 * log.cpue[t - 3]
}
Sigma.max <- 100
m.season <- mean(season.effect[1:N.season])
for (j in 1:N.season) {
season.effect[j] ~ dnorm(0.0, tau[3])
}
for (j in 1:N.sigma) {
tau[j] <- 1 / (sigma[j] * sigma[j])
sigma[j] ~ dunif(0.0, Sigma.max)
}
}