JAGSで予測分布の作成その2

IRTモデルで予測分布や対数尤度の計算などをおこなってみる。
対数密度の記述は以下の頁の質問に少し書かれている。
https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/c728e87b/

The general rule is that if y ~ dfoo(alpha, beta, gamma, ...) then the contribution to the log likelihood of y is given by lik[i] <- logdensity.foo(y, alpha, beta, gamma, ...).

とのことなので，これに従って書いてみよう。(binomial分布はlogdensity.binなので注意）
また，bugsモジュールも読み込んでおく。
(どうでもいいが，T(0,)という関数はサンプリングする範囲を制限するために使うっぽいので，活用すること。)

#
# IRTモデルでの予測分布の作成など。 
#
library("irtoys")
library("jagsUI")
# help(package="jagsUI")
# help(package="irtoys")
#
# 適当な真値
#
Scored2pl

# 適当なデータ
data <- list(x = Scored,
             I = nrow(Scored),
             J = ncol(Scored)
)


# JAGS code 
cat("model{
    # likelihood
    for(i in 1:I){
      for(j in 1:J){
        x[i, j] ~ dbern(prob[i, j])
        logit(prob[i, j]) <- alpha[j] *(theta[i] - delta[j])
        # 予測分布
        x_pred[i, j] ~ dbern(prob[i, j])
        
        # 個々人の個々の項目の尤度
        logLik[i,j]  <- logdensity.bern(x[i,j],prob[i,j])
        #logLik_temp[i,j]  <- logdensity.bern(x[i,j],prob[i,j])
      }
      #logLik[i] <- sum(logLik_temp[i,])

      # prior
      theta[i] ~ dnorm(0.0, 1.0)
    }
    
  # prior
    for(j in 1:J){
      delta[j] ~ dnorm(0 , pr.delta)
      alpha[j] ~ dnorm(0 , pr.alpha)T(0,)
  }
    pr.delta <- pow(100 , -2)
    pr.alpha <- pow(100 , -2)
    }", file="irt_2pl.jag")

fit_jag <- jags(data=data,
                parameters.to.save=c("alpha", "delta","x_pred","logLik","theta"), 
                model.file="irt_2pl.jag",
                modules=c("glm","bugs","basemod"),
                n.chains = 1,
                n.adapt = 0,
                n.iter = 2000,
                n.burnin = 1000,
                n.thin = 2)
fit_jag
Scored2pl$est