機械学習・自然言語処理の勉強メモ

学んだことのメモやまとめ

Stan:階層モデル②

はじめに



前回は階層モデルの基本的なことを書いた。
kento1109.hatenablog.com

今回は実際のデータでもう少し基本を確認する。

過分散



今回は緑本(データ解析のための統計モデリング入門)で紹介されている階層モデルの導入部分を実践する。

データとしては以下を用意する。

x
 [1] 1 8 9 2 4 3 2 9 5 8

緑本では植物の種子数だったが、これは別に何でもよい。
今回は10点満点のテストの生徒の点数などと考える。

これデータの平均と分散を最尤法で求めると以下となる。

n <- length(x)
opt <- optimize(l, interval = c(0, 1), maximum = TRUE)
p <- opt$maximum

p  # mean
[1] 0.5100017

n*p*(1-p)  # var
[1] 2.499

optimizeは最適解を計算してくれる関数。

しかし、実際の平均と分散を求めると以下となる。

mean(x)
[1] 5.1
var(x)
[1] 9.877778

最尤法で求めた分散と大きく乖離していることが分かる。
(これを「過分散」と言う。)

データ全体は二項分布に従っているが、個人差の影響が大きすぎるためこのようなことが起きてしまっている。

階層モデル



では、Stanで個人差を考慮したモデルを作る。
詳しいことはTJOさんのブログが詳しい。
tjo.hatenablog.com

Stanコードは以下の通りである。

data {
	int<lower=0> N;
	int<lower=0> x[N];
}
parameters {
	real beta;
	real r[N];
	real<lower=0> s;
}
transformed parameters {
	real q[N];
	for (i in 1:N)
		q[i] = inv_logit(beta+r[i]);
}
model {
	for (i in 1:N)
		x[i] ~ binomial(N,q[i]);
	for (i in 1:N)
		r[i]~normal(0,s);
}

推定結果は以下のとおりである。

model <- stan_model("hierarchial.stan")
fit <- sampling(model,data=list(N=length(x),x=x))
fit
Inference for Stan model: hierarchial.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta    0.07    0.02 0.63  -1.16  -0.32   0.06   0.45   1.37   889 1.01
r[1]   -1.84    0.02 1.01  -4.05  -2.44  -1.77  -1.15  -0.07  1923 1.00
r[2]    1.16    0.02 0.93  -0.57   0.55   1.12   1.74   3.16  1486 1.00
r[3]    1.74    0.02 1.00   0.01   1.06   1.68   2.34   3.93  2059 1.00
r[4]   -1.29    0.02 0.91  -3.23  -1.83  -1.24  -0.69   0.38  1592 1.00
r[5]   -0.43    0.02 0.82  -2.08  -0.93  -0.41   0.10   1.18  1313 1.00
r[6]   -0.81    0.02 0.86  -2.64  -1.35  -0.77  -0.24   0.80  1356 1.01
r[7]   -1.28    0.02 0.90  -3.24  -1.81  -1.24  -0.67   0.33  1708 1.00
r[8]    1.72    0.02 1.00  -0.01   1.05   1.63   2.32   3.90  1839 1.00
r[9]   -0.06    0.02 0.82  -1.70  -0.59  -0.06   0.47   1.53  1284 1.00
r[10]   1.15    0.02 0.89  -0.52   0.55   1.11   1.71   2.96  1619 1.00
s       1.75    0.01 0.64   0.83   1.30   1.64   2.08   3.28  1835 1.00
q[1]    0.18    0.00 0.11   0.02   0.09   0.16   0.24   0.44  4000 1.00
q[2]    0.75    0.00 0.13   0.48   0.67   0.76   0.85   0.94  4000 1.00
q[3]    0.83    0.00 0.10   0.58   0.77   0.85   0.91   0.98  4000 1.00
q[4]    0.25    0.00 0.12   0.05   0.16   0.24   0.33   0.51  4000 1.00
q[5]    0.42    0.00 0.14   0.17   0.32   0.41   0.52   0.69  4000 1.00
q[6]    0.34    0.00 0.13   0.11   0.24   0.33   0.43   0.63  4000 1.00
q[7]    0.25    0.00 0.12   0.06   0.16   0.24   0.33   0.53  4000 1.00
q[8]    0.83    0.00 0.11   0.57   0.77   0.85   0.91   0.98  4000 1.00
q[9]    0.50    0.00 0.14   0.24   0.40   0.50   0.60   0.77  4000 1.00
q[10]   0.75    0.00 0.12   0.48   0.67   0.76   0.84   0.95  4000 1.00
lp__  -63.37    0.08 2.80 -69.92 -65.00 -63.00 -61.33 -58.96  1340 1.00

これだけではよくわからない。

個人の二項分布の事象の生起確率は以下より、q[i]と分かる。

q[i] = {\rm invlogit}(\beta+r[i])\\
x[i] \sim {\rm binomial}(N,q[i])

これに試行回数Nを掛けると、各個人の確率分布の平均が分かる。
これより、各生徒の二項分布の平均は
1.8, 7.5, 8.3, 2.5, 4.2, 3.4, 2.5, 8.3, 5.0, 7.5」だと分かる。
実際の点数は

x
 [1] 1 8 9 2 4 3 2 9 5 8

だったので、割といい線いっていることが分かる。