SSブログ

JAGSとlooパッケージでWAICを計算してみた [統計]

looパッケージに、WAICを計算する関数waic()があることを最近になってやっと知ったので、つかってみました。

@berobero11さんがStanで計算させた例と、@siero5335さんがStanとlooとで計算させた例とがありますので、それと同じデータで比較してみました。

データを生成します。

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とほぼ同じ値がえられました。


タグ:jags WAIC
nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:日記・雑感

nice! 1

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0