ベイズの定理で二項分布の推定~最尤法との比較まで~ - 医療職からデータサイエンティストへ

前回は、Rにstanを導入して、実際に動かすところまで行いました。今回は、ベイズ推定で二項分布のパラメーターを求めてみます！

medi-data.hatenablog.com

二項分布とは

ベイズの定理を使う前に、さらっと二項分布を復習しておきます。二項分布とは、確率 $\mu$ で表がでるコインをN回投げて、表がm回、裏がN-m回でた場合に、以下のように表すことができる確率分布になります。

$Bin(m|N,\mu) = {}_N C _m \mu^{m}(1-\mu)^{N-m}$

例えば表の裏のでる確率が等しく0.5だった場合に、コインを10回投げて、表が5回でる確率を求める時には、N=10、m=5、 $\mu$ =0.5として計算すると約25%になります。

表の出る回数mを0~10回と変化させて、縦軸を確率としてグラフにするとこんな感じになります。

stan

m=5の時が一番が高く、表が5回でる確率(25%)が一番高いことが分かります！

二項分布についてはこちらをご参照ください。

medi-data.hatenablog.com

つまり、上記の式 $Bin(m|N,\mu)$ の意味は、表がでる確率が $\mu$ のコインをN回投げて、表がm回でる確率という条件付き確率と同じ意味になります。ちなみにこのm回の出やすさを尤度と呼びます。

medi-data.hatenablog.com

二項分布のパラメーター推定

さて、先ほどは確率が0.5と既知の場合の表が出る確率でした。ここで、あるコインを100回投げて表が70回、裏が30回出たと仮定します。

このコインの表のでる確率 $\mu$ はどれぐらいになるでしょうか？

以下では最尤法とベイズ推定の2つの方法を使っていきます。

最尤法

まずは、よく使われる最尤法でパラメーターを求めてみます。先ほどの二項分布の式にN=100、m=70を入れてみます。

$Bin(70|100,\mu)={}_{100} C _{70} \mu^{70} (1-\mu)^{30}$

これを最大化すれば良いので、簡単にするために対数をとると、

$\ln{Bin(70|100,\mu)}=\ln{{}_{100} C _{70} \mu^{70} (1-\mu)^{30}}\\ = \ln {}_{100} C _{70}+70\ln{\mu} + 30\ln(1-\mu)$

となって、これを対数尤度関数と言います。

この対数尤度関数の第一項は定数なので、第二項以降を $\mu$ で偏微分して、傾きを0とすると

$\dfrac{70}{\mu}-\dfrac{30}{1-\mu} = 0$

となるので、対数尤度関数が最大になるように求めた $\mu$ の値はサンプルの平均である0.7になります。

先ほどの対数尤度関数を横軸を $\mu$ 、縦軸を尤度として、グラフ化してみても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%である可能性が高いと言えます。

ベイズの定理

次は最尤法ではなく、ベイズの定理を使って求めてみます。ベイズの定理を使って式で表すと

$p(\mu | N,m) = \dfrac{p(m|N,\mu)p(\mu)}{p(m)}$

となります。
この時、 $p(m)$ は $\mu$ の値と関係ないので、定数として扱うと

$p(\mu | N,m) \propto p(m|N,\mu)p(\mu)$

となります！（ $\propto$ は比例しているという意味です）

この左辺を事後分布と呼び、パラメーターがN=10、m＝70だった場合の $\mu$ の確率分布となるので、一番可能性のある $\mu$ を求めることができます。

次に、この右辺をよくみてみると、 $p(m|N,\mu)$ は,表がでる確率が $\mu$ のコインをN回投げて、表がm回でる確率を表しているので、これはまさに二項分布の尤度になります。

さて、残りの $p(\mu)$ はどうするかというと、これは事前分布とも呼ばれ、ベータ分布に従うと仮定します。詳しい説明は省きますが、ベータ分布は以下のように表します。

$Beta(\mu| a,b) = \dfrac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1}$

すると事後分布は、

$p(\mu| N,m) \propto Bin(m | N,\mu)Beta( \mu|a,b)\\ \propto \mu^{m+a-1}(1-\mu)^{N-m+b-1}$

と表すことができます。

ベータ分布のパラメータa,bが1の時にベータ分布は以下のような[0,1]の一様分布になり（つまり[ $\mu$ はなんでも良い）、上記の式は二項分布の尤度関数と同じ形になります。

準備が整いました。実際に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は求めたいパラメーター（今回は $\mu$ ）を設定します。あらかじめ上限値、下限値が分かっていればそれを指定する必要があります。今回は、確率なので、下限が0、上限が1ですね！

値を設定する際は"int"や"real"など型を宣言する必要があります。今回の場合は、値が離散値なら"int"、連続値なら"real"を宣言します。他にも"vector"や"matrix"などあるようですが、まだまだ勉強不足です。

modelはそれぞれのパラメーターがどんな分布に従うかを記載します。今回は、 $\mu$ は、a=b=1のベータ分布、表が出た回数mは試行回数N、確率 $\mu$ の二項分布とします。

それぞれの行末には、セミコロンをつけるようにしましょう。

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")

f:id:h-wadsworth02:20190104113512p:plain

最尤法で求めた値と、事後分布の平均値はほぼ0.7と同じ値になってますね。

最尤法と何が違うのか

二つのやり方で、コインの表がでる確率を求めてきましたが、ベイズを使っても最尤法を使っても同じ結果なら簡単な最尤法でいいじゃないか！

と思うかもしれませんが、ベイズ推定の良いところは、表がでる確率 $\mu$ が一つの値ではなく、確率分布として求まることです。

例えば先ほどの結果で、このコインは本当は50%の確率で表がでるけど、100回しかコインを投げてないのでたまたまそんな結果になったんじゃないか？

と考えたとします。最尤法では誰がなんと言おうと70％となってしまうので、これでは議論が進みません。

しかし、ベイズでは $\mu$ が確率分布として表現されているため、 $\mu$ が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とするとベータ分布は以下のようになって、 $\mu =0.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