ベータ分布 - 統計学備忘録（R言語のメモ）

完全独習ベイズ統計学入門

作者: 小島寛之
出版社/メーカー: ダイヤモンド社
発売日: 2015/11/20
メディア: 単行本（ソフトカバー）
この商品を含むブログ (6件) を見る

この本を参考にベータ分布を勉強します.
ベータ分布：ベータ関数により導かれる分布. ベイズ推定の際に事前分布として利用します.

$f(x)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)} \ \ \ \ \ (0＜x＜1)$

ベータ関数 $\ \ \ \ B(\alpha,\beta)=\int _0 ^1 x^{\alpha-1}(1-x)^{\beta-1} dt$

期待値 $\ \ \ \frac{\alpha}{\alpha+\beta}$

分散 $\ \ \ \ \frac{\alpha \beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)}$

グラフを描いてみます
ベータを1に固定してαのみ大きくすると、1より大きくなると指数関数となります．

curve(dbeta(x,1,1),from = 0,to=1,ylim=c(0,5))     #黒
par(new=T)
curve(dbeta(x,2,1),from = 0,to=1,ylim=c(0,5),col=2)    #赤
par(new=T)
curve(dbeta(x,3,1),from = 0,to=1,ylim=c(0,5),col=3)    #緑
par(new=T)
curve(dbeta(x,4,1),from = 0,to=1,ylim=c(0,5),col=4)    #青

f:id:yoshida931:20180326140755p:plain:w400

αとβをランダムに

curve(dbeta(x,1,1),from = 0,to=1,ylim=c(0,3))    #黒
par(new=T)
curve(dbeta(x,2,2),from = 0,to=1,ylim=c(0,3),col=2)   #赤
par(new=T)
curve(dbeta(x,3,3),from = 0,to=1,ylim=c(0,3),col=3)    #緑
par(new=T)
curve(dbeta(x,5,3),from = 0,to=1,ylim=c(0,3),col=4)    #青

f:id:yoshida931:20180326141950p:plain:w400

ベータ分布の使用例1

以下は自作（架空）の問題です．
問）ある疾患の患者群で10中で3人が転倒した場合の、転倒確率をベータ分布によるベイズ統計で推測せよ．
10人中3人なので $\frac{3}{10}$ といのが直観的な確率になるのですが、ベイズ推定では次のように考えることができます. まず転倒する確率をx、転倒しない確率を1-ｘとします．転倒する確率(ｘ)+転倒しない確率(1-x)＝1．xは連続無限に存在するので、確率密度となります（0≦x≦1）．ここで、xの事前分布設定のために、xがどの値でも「同様に確からしい」と仮定します．つまりxの事前分布を下図のような一様分布と仮定します．

$y=1\ \ \ \ (0 \leqq x \leqq1)$

f:id:yoshida931:20180326145458p:plain:w300

「転倒した患者3名、転倒しなかった患者7名」という結果になる確率は $P=x^{3} (1-x)^{7}$ となり、この数式は $\alpha=4, \beta=8$ のベータ分布です.
したがって転倒確率は $\ \ \ \frac{\alpha}{\alpha+\beta}=\frac{4}{12}=\frac{1}{3}$ となります．

ベータ分布の使用例2

A君は身体に障害をもっており、ただいま装具装着の練習中です．30秒以内で装着できる成功確率をxとして考えてみます．今日は8回練習して、5回成功しました．A君の装具装着の成功率はどれくらいでしょうか？ベータ分布の使用例1と同じように、xの事前分布設定のために、xがどの値でも「同様に確からしい」と考えて一様分布を仮定します．
$y=1\ \ \ \ (0 \leqq x \leqq1)$
確率は $\ \ \ \frac{\alpha}{\alpha+\beta}=\frac{6}{10}=\frac{3}{5}$ となります．現在のA君の装具装着確率は以下のようになります．

curve(dbeta(x,6,4),from = 0,to=1,ylim=c(0,3))

f:id:yoshida931:20180329150808p:plain:w450

事前分布が同様に確からしいということもありえないので、A君の今日の成功確率を事前分布として確率分布を推定してみます．ここでは正規化定数となるベータ関数を無視して確率分布のみを求めてみます．
P( 現状の成功確率 , 成功確率 ) = P(x) × 事前分布 = $x (x^{5}) (1-x)^{3}=x^{6}(1-x)^{3}$
$\alpha=7, \beta=4$ となり以下のような確率分布（赤）となります．

curve(dbeta(x,6,4),from = 0,to=1,ylim=c(0,3))
par(new=T)
curve(dbeta(x,7,4),from = 0,to=1,ylim=c(0,3),col=2)

f:id:yoshida931:20180329154806p:plain:w450

感想）理学療法士が適切な事前分布を設定して、子どもたちの可能性を探るって重要なことだと思いました．