JAGSでの予測分布の作成

やり方は簡単で，予測分布を未知の変数としてモデルに書けばいい。
Stanではgenerated quantitiyブロックでrng系の関数を使えばいいが，JAGSではそのまま"x_pred ~ dnorm(mu, sigma)"などとするとRでいうところのx_pred <- rnom(1,mu, sigma)のように乱数を生成してくれる。
以下では母相関係数の推定を行ったあとに，2変量のデータの予測分布を作成してみる。

#
# 予測分布の作成
#

library("jagsUI")

#
# 相関係数については，以下の頁のものを参考にした。
# https://kyougokumakoto.blogspot.jp/2016/09/jags_91.html
#
library(MASS)
set.seed(1234) 
mu <- c(0,0) 
Sig <- matrix(c(1,0.75,0.75,1),2,2) 
x <- mvrnorm(100,mu,Sig) 
cor.test(x[,1],x[,2])

data = list(x = x, N = nrow(x)) 

# JAGS code 
cat("model{
    # likelihood
    for(i in 1:N){
    x[i,1:2] ~ dmnorm(mu[], prec[ , ]) 
    
    # 対数尤度
    logLik[i] <- logdensity.mnorm(x[i,1:2], mu, prec) 
    } 
    
    # 予測分布
    x_pred[1:2] ~ dmnorm(mu[], prec[ , ]) 

    # corresponding the precision matrix
    prec[1:2, 1:2] <- inverse(cov[,]) 
    
    # constructing the covariance matrix
    cov[1,1] <- sigma[1]^2
    cov[1,2] <- sigma[1]*sigma[2]*rho
    cov[2,1] <- sigma[1]*sigma[2]*rho
    cov[2,2] <- sigma[2]^2 
    
    #priors
    rho ~ dunif(-1,1) 
    mu[1] ~ dnorm(0, 1.0e-4) 
    mu[2] ~ dnorm(0, 1.0e-4)
    sigma[1] ~ dt(0, 1, 1)T(0,)
    sigma[2] ~ dt(0, 1, 1)T(0,) 
    
    }", file="cor_pred.jag")

fit_jag <- jags(data=data,
                parameters.to.save=c("mu", "sigma","rho","x_pred","logLik"), 
                model.file="cor_pred.jag",
                n.chains = 1,
                n.adapt = 100,
                n.iter = 2000,
                n.burnin = 1000,
                n.thin = 2)

# 予測分布の相関
cor(fit_jag$samples[[1]][,c("x_pred[1]","x_pred[2]")])
# 予測分布のヒストグラムの重ね合わせ
hist(fit_jag$samples[[1]][,c("x_pred[1]")],freq = F,col=rgb(1,0,1,alpha=0.5))
hist(x[,1],freq = F,add=T,col=rgb(1,0,0,alpha=0.5))