前回は、Rにstanを導入して、実際に動かすところまで行いました。今回は、ベイズ推定で二項分布のパラメーターを求めてみます!
二項分布とは
ベイズの定理を使う前に、さらっと二項分布を復習しておきます。 二項分布とは、確率で表がでるコインをN回投げて、表がm回、裏がN-m回でた場合に、以下のように表すことができる確率分布になります。
例えば表の裏のでる確率が等しく0.5だった場合に、コインを10回投げて、表が5回でる確率を求める時には、N=10、m=5、=0.5として計算すると約25%になります。
表の出る回数mを0~10回と変化させて、縦軸を確率としてグラフにするとこんな感じになります。
m=5の時が一番が高く、表が5回でる確率(25%)が一番高いことが分かります!
二項分布についてはこちらをご参照ください。
つまり、上記の式の意味は、表がでる確率がのコインをN回投げて、表がm回でる確率という条件付き確率と同じ意味になります。 ちなみにこのm回の出やすさを尤度と呼びます。
二項分布のパラメーター推定
さて、先ほどは確率が0.5と既知の場合の表が出る確率でした。ここで、あるコインを100回投げて表が70回、裏が30回出たと仮定します。
このコインの表のでる確率はどれぐらいになるでしょうか?
以下では最尤法とベイズ推定の2つの方法を使っていきます。
最尤法
まずは、よく使われる最尤法でパラメーターを求めてみます。先ほどの二項分布の式にN=100、m=70を入れてみます。
これを最大化すれば良いので、簡単にするために対数をとると、
となって、これを対数尤度関数と言います。
この対数尤度関数の第一項は定数なので、第二項以降をで偏微分して、傾きを0とすると
となるので、対数尤度関数が最大になるように求めたの値はサンプルの平均である0.7になります。
先ほどの対数尤度関数を横軸を、縦軸を尤度として、グラフ化してみても0.7のところが最大になっていることが分かります。
>#対数尤度関数の定義 > log_like <- function(mu){ 70*log(mu) + 30*log(1-mu) } > mu <- seq(0,1,0.1)#確率を0~1まで0.1づつ変化させる > plot(mu,log_like(mu),type = "l" ,xlab = "",ylab ="",main="" ,col = "red",lwd =2,xaxt = "n") > abline(v = 0.7) > axis(1,mu,mu) >
つまりこのコインは表のでる確率が70%である可能性が高いと言えます。
ベイズの定理
次は最尤法ではなく、ベイズの定理を使って求めてみます。ベイズの定理を使って式で表すと
となります。
この時、はの値と関係ないので、定数として扱うと
となります!(は比例しているという意味です)
この左辺を事後分布と呼び、パラメーターがN=10、m=70だった場合のの確率分布となるので、一番可能性のあるを求めることができます。
次に、この右辺をよくみてみると、は,表がでる確率がのコインをN回投げて、表がm回でる確率を表しているので、これはまさに二項分布の尤度になります。
さて、残りのはどうするかというと、これは事前分布とも呼ばれ、ベータ分布に従うと仮定します。詳しい説明は省きますが、ベータ分布は以下のように表します。
すると事後分布は、
と表すことができます。
ベータ分布のパラメータa,bが1の時にベータ分布は以下のような[0,1]の一様分布になり(つまり[はなんでも良い)、上記の式は二項分布の尤度関数と同じ形になります。
準備が整いました。実際にstanを使ってみます。まずは、stanファイルを新規で作成して、以下を"binominal.stan"として保存します。
data { int N;//試行回数 int m ;//表がでた回数 int a ;//ベータ分布のパラメーター int b ; } parameters { real<lower=0,upper=1> mu; } model { mu~beta(a,b); m~binomial(N,mu); }
この辺りの記載方法は少し難しいですが、慣れるしかなさそうです。
dataはRスクリプトからデータを渡す時に必要な項目です。
parametersは求めたいパラメーター(今回は)を設定します。あらかじめ上限値、下限値が分かっていればそれを指定する必要があります。今回は、確率なので、下限が0、上限が1ですね!
値を設定する際は"int"や"real"など型を宣言する必要があります。今回の場合は、値が離散値なら"int"、連続値なら"real"を宣言します。他にも"vector"や"matrix"などあるようですが、まだまだ勉強不足です。
modelはそれぞれのパラメーターがどんな分布に従うかを記載します。今回は、は、a=b=1のベータ分布、表が出た回数mは試行回数N、確率の二項分布とします。
それぞれの行末には、セミコロンをつけるようにしましょう。
Rのスクリプトファイルには、以下のように記載します。
>library(rstan) > #stanファイルを読み込む > stanmodel <- stan_model("binominal.stan") DIAGNOSTIC(S) FROM PARSER: Info (non-fatal): Comments beginning with # are deprecated. Please use // in place of # for line comments. Info (non-fatal): Comments beginning with # are deprecated. Please use // in place of # for line comments. Info (non-fatal): Comments beginning with # are deprecated. Please use // in place of # for line comments. >#サンプリングする > paralist <- list(N=100,m=70,a=1,b=1) >fit <- sampling(stanmodel,data = paralist)#パラメーターはリストで渡す >> fit Inference for Stan model: binominal. 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 mu 0.70 0.00 0.05 0.60 0.67 0.70 0.73 0.78 1575 1 lp__ -63.17 0.02 0.76 -65.34 -63.34 -62.89 -62.70 -62.64 1503 1 Samples were drawn using NUTS(diag_e) at Thu Jan 3 17:24:42 2019. 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).
おお!初Stanでしたがうまくいったようです!結果を抜き出すには、rstan::extraを使います。
> result <- rstan::extract(fit) > mean(result$mu) [1] 0.6965801 > hist(result$mu,freq = F,main ="" + ,xlab ="",ylab ="") >lines(density(result$mu),col ="red")
最尤法で求めた値と、事後分布の平均値はほぼ0.7と同じ値になってますね。
最尤法と何が違うのか
二つのやり方で、コインの表がでる確率を求めてきましたが、ベイズを使っても最尤法を使っても同じ結果なら簡単な最尤法でいいじゃないか!
と思うかもしれませんが、ベイズ推定の良いところは、表がでる確率が一つの値ではなく、確率分布として求まることです。
例えば先ほどの結果で、このコインは本当は50%の確率で表がでるけど、100回しかコインを投げてないのでたまたまそんな結果になったんじゃないか?
と考えたとします。最尤法では誰がなんと言おうと70%となってしまうので、これでは議論が進みません。
しかし、ベイズではが確率分布として表現されているため、が0.5になる確率を求めることができます。
Rで実行してみると
>#muが0.5以上、0.6未満となる確率 > sum((result$mu>=0.5) & (result$mu<0.6))/length(result$mu) [1] 0.02175
となりおよそ2%の確率で、表がでる確率が50%であることが分かりました。
さらに、事前分布のベータ分布を一様分布としましたが、a=b=5とするとベータ分布は以下のようになって、の時が一番大きくなります。
これを事前分布とすると、コインは50%の確率で表が出るはず!という前提をモデルに組み込むことができます。
>#事前分布にa=b=5のベータ分布を設定 > paralist <- list(N=100,m=70,a =5,b=5) > fit2 <- sampling(stanmodel,data = paralist) > result2 <- rstan::extract(fit2) > plot(density(result2$mu)$x,density(result2$mu)$y + ,xlab ="",ylab ="",type ="l",col ="red") > lines(density(result$mu)) > > sum((result2$mu>=0.5) & (result2$mu<0.6))/length(result2$mu) [1] 0.0335
グラフをみてわかるように、今回(赤線)の方が全体的に分布が左によっています。50%になる確率も3.4%と若干あがりました。
今回の例のように、ベイズの定理を使うと事後確率であるパラメーターの分布から多くの情報を得ることができたり、事前分布を変えることで、前提条件をモデルに組み込むこともできます!
これが最尤法にはないメリットですね。
まとめ
今回はstanの実践練習ということで、ベイズ推定を使って簡単なパラメーター推定を行いました。
しかし、stanのすごいところ最尤法が使えないようなモデリングも行うことができるところにあります!
今後は、階層ベイズモデルなどにも挑戦していきます。
※本記事は筆者が個人的に学んだことをまとめた記事なります。数学の記法や詳細な理論、筆者の勘違い等で誤りがあった際はご指摘頂けると幸いです。
参考
ベイズモデリングをする方なら必読書と言っても過言ではありません。これほど分かりやすく学べる参考書はないでしょう。
StanとRでベイズ統計モデリング (Wonderful R)
- 作者: 松浦健太郎,石田基広
- 出版社/メーカー: 共立出版
- 発売日: 2016/10/25
- メディア: 単行本
- この商品を含むブログ (10件) を見る