JAGSとlooパッケージでWAICを計算してみた [統計]
looパッケージに、WAICを計算する関数waic()があることを最近になってやっと知ったので、つかってみました。
@berobero11さんがStanで計算させた例と、@siero5335さんがStanとlooとで計算させた例とがありますので、それと同じデータで比較してみました。
- @berobero11さん:[Stan] WAICとWBICを事後分布から計算してみる
- @siero5335さん:loo package動かしてみた: WAIC比較
データを生成します。
N <- 100 a_true <- 0.4 mean1 <- 0 mean2 <- 3 sd1 <- 1 sd2 <- 1 set.seed(1) Y <- c(rnorm((1 - a_true) * N, mean1, sd1), rnorm(a_true * N, mean2, sd2))
Model 1は混合分布モデルなのですが、JAGSにはincrement_log_prob()がないので、潜在離散パラメータ(z)を使用してモデリングしています。ただし、対数尤度の計算の部分はStanと同じにしています。
## model1a.bug.txt var # Data N, # Sample size Y[N], # Observed data # Parameters a, # Mixing rate mu, # Mean for the 2nd group z[N], # Latent parameter loglikelihood[N]; model { # Priors a ~ dunif(0, 1); mu ~ dunif(-50, 50); # Likelihood for (i in 1:N) { z[i] ~ dbern(a); Y[i] ~ dnorm(z[i] * mu, 1); } # Generated quantities for (i in 1:N) { loglikelihood[i] <- log((1 - a) * dnorm(Y[i], 0, 1) + a * dnorm(Y[i], mu, 1)); } }
Model 2は単純に正規分布にあてはめています。
## model2a.bug.txt var # Data N, # Sample size Y[N], # Observed data # Parameters mu, # Mean s, # SD loglikelihood[N]; model { # Priors mu ~ dnorm(0, 1.0e-4); s ~ dunif(0, 1000); # Likelihood for (i in 1:N) { Y[i] ~ dnorm(mu, 1 / (s * s)); } # Generated quantities for (i in 1:N) { loglikelihood[i] <- log(dnorm(Y[i], mu, 1/ (s * s))); } }
rjagsからJAGSを実行します。
library(rjags) ## Model 1 model1 <- jags.model("model1a.bug.txt", data = list(N = N, Y = Y), n.chains = 4, n.adapt = 500) update(model1, 500) fit1 <- coda.samples(model1, n.iter = 10000, thin = 10, variable.names = c("a", "mu", "loglikelihood")) ## Model 2 model2 <- jags.model("model2a.bug.txt", data = list(N = N, Y = Y), n.chains = 4, n.adapt = 500) update(model2, 500) fit2 <- coda.samples(model2, n.iter = 1000, thin = 1, variable.names = c("mu", "s", "loglikelihood"))
計算された対数尤度を行列の形に整形して、looパッケージのwaic()関数でWAICを計算します。
library(loo) loglik1 <- sapply(1:N, function(i) unlist(fit1[, paste("loglikelihood[", i, "]", sep = "")])) waic(loglik1) loglik2 <- sapply(1:N, function(i) unlist(fit2[, paste("loglikelihood[", i, "]", sep = "")])) waic(loglik2)
結果です。
> library(loo) This is loo version 0.1.2 > > loglik1 <- sapply(1:N, function(i) + unlist(fit1[, paste("loglikelihood[", i, "]", sep = "")])) > waic(loglik1) Computed from 4000 by 100 log-likelihood matrix Estimate SE Estimate SE elpd_waic -191.4 5.5 p_waic 2.0 0.3 waic 382.7 11.0 > > loglik2 <- sapply(1:N, function(i) + unlist(fit2[, paste("loglikelihood[", i, "]", sep = "")])) > waic(loglik2) Computed from 4000 by 100 log-likelihood matrix Estimate SE elpd_waic -198.1 5.6 p_waic 1.6 0.2 waic 396.1 11.1
Model 1,2ともに、Stanとほぼ同じ値がえられました。
コメント 0