So-net無料ブログ作成

事前分布とWAIC [統計]

2項分布などの確率パラメータに、確率スケールで一様事前分布をつけるか、ロジットスケールで一様事前分布をつけるかで、事後分布が変わるということが指摘されていますが1、WAICでどのくらいの差がでるのか実験してみました。

1 Kéry & Schaub “Bayesian Population Analysis using WinBUGS” 2章とか。

今回はRStanをつかいました。

モデルは以下の2つです。ともにStanの暗黙の事前分布を使用していますが、1つ目のモデルは確率スケールでの一様事前分布、2つ目のモデルはロジットスケールでの非正則一様事前分布となっています。

## binomial1.stan
data {
  int<lower=0> N;
  int<lower=0> M;
  int<lower=0,upper=M> y[N];
}

parameters {
  real<lower=0,upper=1> p;
}

model {
  y ~ binomial(M, p);
}

generated quantities {
  vector[N] log_lik;

  for (i in 1:N)
    log_lik[i] <- binomial_log(y[i], M, p);
}
## binomial2.stan
data {
  int<lower=0> N;
  int<lower=0> M;
  int<lower=0,upper=M> y[N];
}

parameters {
  real logit_p;
}

model {
  y ~ binomial_logit(M, logit_p);
}

generated quantities {
  real p;
  vector[N] log_lik;

  p <- inv_logit(logit_p);
  for (i in 1:N)
    log_lik[i] <- binomial_logit_log(y[i], M, logit_p);
}

Rのコードは以下のようなものをつかいました。小サンプルと大サンプルの2種類のデータを用意しています。

## Use RStan and loo
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(loo)

## Model files
model_file1 <- "binomial1.stan"
model_file2 <- "binomial2.stan"

## Data 1: small sample
N1 <- 16
M <- 10
p <- 0.3
set.seed(3)
y1 <- rbinom(N1, M, p)

fit11 <- stan(model_file1, data = list(N = N1, M = M, y = y1),
              iter = 2000, warmup = 1000, chains = 4,
              open_progress = FALSE)
fit12 <- stan(model_file2, data = list(N = N1, M = M, y = y1),
              iter = 2000, warmup = 1000, chains = 4,
              open_progress = FALSE)
## Summary
print(fit11, pars = "p", digits = 3)
print(fit12, pars = "p", digits = 3)

## WAIC
loglik11 <- extract_log_lik(fit11)
loglik12 <- extract_log_lik(fit12)
waic(loglik11)
waic(loglik12)

## Data 2: large sample
N2 <- 400
M <- 10
p <- 0.3
set.seed(31)
y2 <- rbinom(N2, M, p)

fit21 <- stan(model_file1, data = list(N = N2, M = M, y = y2),
              iter = 2000, warmup = 1000, chains = 4,
              open_progress = FALSE)
fit22 <- stan(model_file2, data = list(N = N2, M = M, y = y2),
              iter = 2000, warmup = 1000, chains = 4,
              open_progress = FALSE)

## Summary
print(fit21, pars = "p", digits = 3)
print(fit22, pars = "p", digits = 3)

## WAIC
loglik21 <- extract_log_lik(fit21)
loglik22 <- extract_log_lik(fit22)
waic(loglik21)
waic(loglik22)

結果です。

> waic(loglik11)
Computed from 4000 by 16 log-likelihood matrix

          Estimate  SE
elpd_waic    -24.6 1.3
p_waic         0.4 0.2
waic          49.2 2.6
> waic(loglik12)
Computed from 4000 by 16 log-likelihood matrix

          Estimate  SE
elpd_waic    -24.6 1.3
p_waic         0.4 0.2
waic          49.1 2.6
> waic(loglik21)
Computed from 4000 by 400 log-likelihood matrix

          Estimate   SE
elpd_waic   -714.4 12.7
p_waic         1.0  0.1
waic        1428.8 25.4
> waic(loglik22)
Computed from 4000 by 400 log-likelihood matrix

          Estimate   SE
elpd_waic   -714.4 12.7
p_waic         1.0  0.1
waic        1428.8 25.4

この単純な例では、WAICにはほとんど差はありませんでした。

これだけではあまりおもしろくないので、情報事前分布をつけたときどうやるかもやってみました。logit_pにNormal(5, 1)という事前分布を設定しています。

## binomial3.stan
data {
  int<lower=0> N;
  int<lower=0> M;
  int<lower=0,upper=M> y[N];
}

parameters {
  real logit_p;
}

model {
  logit_p ~ normal(5, 1);
  y ~ binomial_logit(M, logit_p);
}

generated quantities {
  real p;
  vector[N] log_lik;

  p <- inv_logit(logit_p);
  for (i in 1:N)
    log_lik[i] <- binomial_logit_log(y[i], M, logit_p);
}

結果です。

> waic(loglik13)
Computed from 4000 by 16 log-likelihood matrix

          Estimate  SE
elpd_waic    -25.0 1.3
p_waic         0.4 0.2
waic          50.0 2.6
> waic(loglik23)
Computed from 4000 by 400 log-likelihood matrix

          Estimate   SE
elpd_waic   -714.4 12.7
p_waic         1.0  0.1
waic        1428.8 25.4

小サンプルではWAICがややおおきくなりましたが、大サンプルではほとんど変化ありませんでした。サンプルサイズがおおきくなると事前分布の影響はちいさくなりますので、まあ当然の結果です。


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

nice! 1

コメント 0

コメントを書く

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

Facebook コメント

トラックバック 0