理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

二項分布の最尤推定法 

  母集団Aの成功率0.7をパラメータとして、この母集団のサンプルから考えていきます.成功率0.7、試行回数10回の二項分布の乱数を1000個生成して母集団Aとします.

y <- rbinom(1000,10,0.7)
head(y)   #最初の6データを確認してみます
[1] 6 7 6 4 9 8  #成功率0.7なので7に近い値になっています
hist(y, xlim=c(0,10))    #ヒストグラムを作成してみます

母集団Aの成功率は0.7なので、10回試行で7回成功したデータが最も多く生成されることになります.
f:id:yoshida931:20180110135858p:plain:w400

次に、試行回数10回で成功率0.3(黒)、0.5(赤)で、成功回数0~10回の確立密度のグラフを描いてみます

x <- rep(0:10)
prob <- dbinom(x,10,0.3)
plot(prob,type="b",xaxt="n")
axis(side=1,at=1:11,labels=x)
prob2 <- dbinom(x,10,0.5)
lines(prob2,type="b",xaxt="n",col="2")    #赤色

f:id:yoshida931:20180110140135p:plain:w400
  やはりどちらも母集団Aの分布とは異なります.本来であればパラメータは未知な数値です.どうのようにしてサンプルから推定するのでしょうか?.実際に10回中7回が成功したら, 成功確立は 7 ÷ 10 = 0.7となり, 真のパラメートも0.7になりそうですが…


最尤原理
  最尤原理は、「現在起きている事象の起きる確率」が最大値となるとき、母数の推定値として「最も尤もらしい」とする考え方です.

  例)成功率0.7の母集団から10回試行した場合の成功回数のサンプルを5つ用意して、「現在起きている事象」とします.

#母集団から抽出したサンプルを仮定します
set.seed(5)
rb <- rbinom(5,10,0.7)
[1] 8 6 5 8 9                    #現在起きている事象=成功回数

  試行回数は10回なので、それぞれの事象の起きる確率Y_iは、「0.8, 0.6, 0.5, 0.8, 0.9」 となります.
  Y_1=0.8,  Y_2=0.6,  Y_3=0.5,  Y_4=0.8,  Y_5=0.9

母集団の確率密度関数 f\theta(y)=_nC_y{\theta}^y{(1-\theta)}^{n-y}

#母集団から抽出したサンプルの確率密度
#これは成功率を固定して、成功回数が変化しているので確率です
choose(10,8)*0.7^8*0.3^2
choose(10,6)*0.7^6*0.3^4
choose(10,5)*0.7^5*0.3^5
choose(10,8)*0.7^8*0.3^2
choose(10,9)*0.7^9*0.3^1

#Rの関数を使用
dbinom(rb,10,0.7)

[1] 0.2334744 0.2001209 0.1029193 0.2334744 0.1210608


母集団から抽出したサンプルの同時確率密度
サンプル(測定値)「8, 6, 5, 8, 9 」は、互いに独立しています.したがって、「8, 6, 5, 8, 9」が揃う確率は,

f\theta(y_1,y_2,...,y_n)= \prod_{i=0}^n f\theta (y_i)

という同時確率になり、最尤推定値を求めるときに必要になります.

#Rを使用して、母集団から抽出したサンプルの同時確率密度を求めてみます
prod(dbinom(rb,10,0.7))
[1] 0.0001359164
#母集団の成功率0.7が前提だったので算出可能となります


ここからが本題です

  成功回数Y_ia_1, a_2, … ,a_iが観察されたとき、その同時確率密度を母数(成功率)\thetaの関数とみて、尤度関数といいます. 尤度関数は、成功回数(y1=8,  y2=6 , y3=5 , y4=8 , y5=9 )を固定して、成功率\thetaが変化する関数ということになります.つまり、データを定数として母数を変数とした、「観測されたデータが得られる確率」を「未知のパラメータの関数」とみた関数です.

尤度関数 L(θ)=f\theta(a_1,a_2,...,a_n)= \prod_{i=0}^n f\theta (a_i)

# 測定値「8,6,5,8,9」, 母数θは未定
L(θ) = choose(10,8)*θ^8*(1-θ)^2 ×
        choose(10,6)*θ^6*(1-θ)^4 ×
        choose(10,5)*θ^5*(1-θ)^5 ×
        choose(10,8)*θ^8*(1-θ)^2 ×
        choose(10,9)*θ^9*(1-θ)^1



最尤推定
  上述の例題の\thetaは0.7だったのですが、母集団の母数が分かっていることはほとんどありません.本来であれば不明です.したがって、成功率\thetaが不明の場合は、尤度関数から\thetaを求めます.

尤度関数 L(θ)= \prod_{i=0}^n f\theta (a_i)

  \thetaを尤度関数から推定する場合に、尤度関数を最大にする\thetaa_1, a_2, … ,a_iが観測されたときの最尤推定値といいます. したがって、尤度関数を微分して L'(θ)=0 となるθを求めます.
  とはいっても、尤度関数から直接計算するのは、なんとなく大変そうなので、対数をとって考えていきます.

二項分布の対数尤度関数
上述の例題では試行回数n=10回、成功回数a_i   (i = 1,2,3,4,5    )

f\theta(ai) =_{10}C_a{_i}{\theta}^{ai}{(1-\theta)}^{10-ai}

=log[\frac{n!}{x!(n−x)!}θ^{ai}(1−θ)^{n−x}]

=log(n!)−log(x!)−log(n−x)!+ailogθ+(10-ai)log(1−θ)

対数尤度関数を微分
=log(n!)−log(x!)−log(n−x)!をθで微分すると0

したがって f' \theta(ai) = \sum(\frac{ai}{\theta}-\frac{10-ai}{1−θ})

ai=(8,6,5,8,9)を挿入

f' \theta(ai) =\frac{8}{\theta}-\frac{10-8}{1−θ}+\frac{6}{\theta}-\frac{10-6}{1−θ}+\frac{5}{\theta}-\frac{10-5}{1−θ}+\frac{8}{\theta}-\frac{10-8}{1−θ}+\frac{9}{\theta}-\frac{10-9}{1−θ}=0
  = \frac{36-50\theta}{\theta(1-\theta)}=0
\theta=0.72
パラメータは0.72と推定することができました.

#Rを使って確認してみます
LL <- function(p) {
  prod(dbinom(x,10,p))           # 同時密度分布
}
optimize(f=LL, c(0, 1), maximum=TRUE)      #optimizeでパラメータの最適化(この場合は最大値を求める)
$maximum
[1] 0.7200026

参考)合成関数の導関数
y = f(u)  ,  u = g(x) のとき  y = f( g(x) ) となる y = f( g(x) ) 導関数は、
\frac{dy}{dx}=\frac{dy}{du}*\frac{du}{dx}


参考web
norimune.net
www.yasuhisay.info

参考文献
柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010
東京大学教養学部統計学教室 (編集); 統計学入門 (基礎統計学), 東京大学出版会, 1991