コイン投げの事後分布を発生させてみる

データ

我々が入手できるデータは、3回のコイン投げで(表、表、裏)が出たという事実だけです。
これを、次のようなリストとして保存します。

d01<-list(N=3,y=c(1,1,0))

モデルを定義する

このコインの表が出る確率$p$の事後分布を発生させるために、Stanのモデルを記述します。

RMarkdownを使っているなら(ていうか、使いましょう!)、以下を


```{stan,output.var="(スクリプト名)"}

   (Stan のスクリプト)

```

で囲って記述します。


```{stan,output.var="coin01"}
data{
  int N;
  int y[N];
}

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

model{
  for(i in 1:N){
    y[i]~binomial(1,p);
  }
}
```

解説します。

data

最初のdata{...}のところで、読み込むデータがどんなものか、Stanに教えてあげます。
intというのは、「整数だよ」って言っています。
int N;で、「Nは整数だよ」と言っているわけです。
最後の;は記述の最後に付けることになっています。
int y[N];では、「yは整数」だけどこれはN個のベクトル(数の集まり)だということを表しています。

parameters

parameters{...}のところでは、事後分布を推定するパラメータを定義します。正確には、事後分布を推定しなくても、あとのmodelの記述で出てくる未知数のことです。今回は表が出る確率$p$ですね。

realというのは実数(連続変数のこと)です。
その後に<lower=0>と書くと0以上<upper=1>と書くと1以下に分布を制限してくれます。ここでは両方必要なので<lower=0,upper=1>と書いてあります。

model

model{...}のところに尤度を決める確率モデルを定義します。
y[1]~binomial(1,p);で、yの1番目のデータで、1回投げでy[1]が出るのは次の確率であると記述しています。

$$\Pr(y_1=1|p)=p(1-p)$$

データyの要素はN個あるので、for(i in 1:N){...}を使って、y[2]y[3]についても定義しています。
Rだとベクトルを一般に扱ったりできるのですが、Stanでは一個一個定義しないといけません。

また、モデルで使った変数は、必ずdata{...}parameters{...}のところで定義しておかなければなりません。

ちなみに・・・・binomialという関数について説明しておきます。
これは二項分布のことで、コイン投げで、$N$回投げて$n$回表が出る確率を定義します。

$$\Pr(n|N,p)={}_NC_n p^n(1-p)^{N-p}$$

今回は、これを$N=1$(1回だけの試行)、$n$は0か1かがyで与えられるものとして、わざわざ3回分の確率モデルを定義しています。

これはこれまでの説明に沿って解説するためで、もうちょっと簡単にするならば、

model{
     n~binomial(N,p);
}

でもよいわけです。そうなると、data{...}の箇所は、

data{
        int N;
        int n;
}

となります。与えるデータは3回投げて2回が当たりという情報だけで以下となります。
d01<-list(N=3,n=2)
Stanにはさまざまな確率モデルが用意されていて、
http://mc-stan.org/users/documentation/
のDocumentationのVIIIとIXで確認することができます。

ちょっと長くなりましたが、モデルを記述したら、RStudioのsource画面の右上のCheckをクリックして、モデルの記述に間違いがないかチェックします。ただし、ここでチェックできるのは文法的な間違い(;が抜けているとか)だけで、モデルそのものが間違っているかどうかはチェックできません。

実行

ライブラリrstanを読み込みます。すでに読み込んでいればここは不要です。

library(rstan)

いよいよサンプリングをして、事後分布を発生させます。

(fit01<-sampling(object = coin01,
                 data=d01))
| 
| SAMPLING FOR MODEL '2aac5ce27e67f00e55e2cf689479a4cc' NOW (CHAIN 1).
| Chain 1: 
| Chain 1: Gradient evaluation took 0 seconds
| Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
| Chain 1: Adjust your expectations accordingly!
| Chain 1: 
| Chain 1: 
| Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
| Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
| Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
| Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
(中略)
| Inference for Stan model: 2aac5ce27e67f00e55e2cf689479a4cc.
| 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
| p     0.60    0.01 0.2  0.19  0.46  0.61  0.76  0.93  1382 1.01
| lp__ -3.92    0.02 0.8 -6.08 -4.11 -3.62 -3.42 -3.37  1426 1.00
| 
| Samples were drawn using NUTS(diag_e) at Tue Sep 01 22:26:14 2020.
| For each parameter, n_eff is a crude measure of effective sample size,
| and Rhat is the potential scale reduction factor on split chains (at 
| convergence, Rhat=1).

出力の最後の方にサマリーが出ています。上段のpのところが、$p$の事後分布です。

収束の確認

結果が出たからと、すぐに喜んではいけません。ちゃんと分布が収束しているかどうか確認しなければなりません。

事後分布の推定は、定義したモデルにうまく合致するような乱数を発生させて行います(MCMC, マルコフチェーン・モンテカルロと言う方法です)。

発生させた乱数が、ちゃんとある分布に落ち着いたか(定常分布と言います)収束を確認します。

トレースラインの確認

その方法のひとつが、以下のように乱数をを発生させた軌跡をプロットしてみることです。

トレースラインの出力
stan_trace(fit01)

このようにグチャグチャに混ざっているのはうまくできた結果です。
うまくできていないと以下のようになります。前のプロットに次のプロットがひきずられるランダムウォークのような格好となっています。

自己相関の確認

以前のプロットとの相関(自己相関と言います)を確認すると、うまくいった方は、こんな感じです。

自己相関を確認する
stan_ac(fit01,pars="p")

うまくいかなかったのはこんな感じです。

「Lag」というのは何期前かということで、うまくいった方は、4期ぐらい以前にはほとんど相関がないが、うまくいかなかった方は20期前でも結構相関がみられます。

Chain間の差異の確認

プロットの図にChain
1、2、3、4とあるのは、乱数発生の最初に与えた値(初期値)によって、収束する分布が違わないかを確認するため、別々の初期値を与えてみて同じ分布に収束するかを確認します。

Chainごとの分布の違いの確認
stan_dens(fit01,pars="p",separate_chains = T)

ほぼ同じ分布となっていますね。これがうまくいかないと次のようにバラバラの分布となります。

Warmup期間の確認

適当に設定した初期値に引きずられれるので、その期間は捨てます。この期間をBarn-in(焼き入れ)またはWarmup期間と言います。

Warmup期間が十分か確認する
stan_trace(fit01,inc_warmup=TRUE)

これは単純で値域が限られている確率モデルなのですぐに定常分布に収束していますが、通常はいろいろ遠回りることもあります。

図だけでなく、何か判定する指標が欲しかったら、Gelman & Rubin(1992)の$\hat{R}$指標があります。先に示したStanの出力の最後の列Rhatがこれにあたります。この値が 1.1を下回ることが定常分布に収束したことの目安とされています。

事後分布の出力

以上で定常分布への収束を確認したら、事後分布を出力します。収束していなくても事後分布は出力できますが、この分布に意味がありません。定常分布への収束が確認された場合のみにおいて、この作業を行います。

事後分布(ヒストグラム)の出力
stan_hist(fit01,pars="p")
| stat_bin() using bins = 30. Pick better value with binwidth.

事後分布(密度分布)の出力
stan_dens(fit01,pars="p")

サマリーの出力
summary(fit01)$summary
|            mean     se_mean        sd      2.5%        25%        50%
| p     0.5989586 0.005395938 0.2006227  0.188703  0.4569711  0.6113394
| lp__ -3.9218842 0.021149290 0.7985132 -6.078782 -4.1077644 -3.6162753
|             75%      97.5%    n_eff     Rhat
| p     0.7568167  0.9310394 1382.376 1.007396
| lp__ -3.4200723 -3.3655562 1425.518 1.000632

左から、
mean: 平均、
se_mean: 平均の標準誤差、
sd: 標準偏差、
xx%: xx%分位点、
n_eff: 有効プロット数、
Rhat: $\hat{R}$

平均も中央値も0.6ぐらい、2/3であることがわかります。

事前分布を指定する

結果が尤度の分布$\Pr(y|p)=p^2(1-p)$とほとんど一致していることに気づいたでしょうか?それもそのはずで、Stanは、事前分布に何も指定しなかったら、事前分布として一様分布を使います。一様分布の確率密度は、

$$\frac{d}{dp}\Pr(\le p)=1$$

なので、事後分布は

$$\frac{d}{dp}\Pr(\le p|y)=p^2(1-p)\times 1$$

であり、尤度関数そのものです。

事前分布として、平均0.5、標準偏差0.2の正規分布を指定してみましょう。

stanでこの事前分布を指定するのは、model{...}の部分に以下を追加します。

p~normal(0.5,0.2);

すなわち、coin01を以下のように書き換えて、coin02とします。


```{stan,output.var="coin02"}
data{
  int N;
  int y[N];
}

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

model{
  p~normal(0.5,0.2);
  for(i in 1:N){
    y[i]~binomial(1,p);
  }
}
```

そこでもう一度MCMCを実行します。

fit02<-sampling(object=coin02,d01)

結果は、

summary(fit02)$summary
|            mean     se_mean        sd       2.5%        25%        50%
| p     0.5502292 0.004244351 0.1560872  0.2409435  0.4415744  0.5555066
| lp__ -4.0138739 0.025204203 0.8273223 -6.3028777 -4.2142791 -3.6973422
|             75%      97.5%    n_eff     Rhat
| p     0.6607516  0.8455139 1352.421 1.000299
| lp__ -3.4843863 -3.4217267 1077.466 1.002612

となり、密度分布は以下となります。

事前分布が一様分布のときよりも少し左に寄っていますね。

事後分布を確率密度関数で表すと、正規分布の確率密度関数が

$$f(x)=\frac{1}{\sqrt{2pi \sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

なので、$x$を$p$で置き換えて、

$$\frac{d}{dp}\Pr(le p)\propto p^2(1-p)\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(p-\mu)^2}{2\sigma^2}\right)$$

となります。ただし、$\mu=0.5$、$\sigma=0.2$です。
グラフに表すとこんな感じです。

Stanで推定した事後分布とほぼ一致してますね。
こんな数式で事後分布がわかるのなら、わざわざStan使ってMCMCやる必要ないじゃないかと言われそうですが、今回はものすごく単純な確率モデルでしたので、こんなに簡単に事後分布が出せただけです。実践的にはこんなもっと複雑で、指定したい事後分布の分布形を解析的に解くのが難しい場合も多いのです。