能力パラメータの推定

1 最尤推定法.

尤度関数.

\(\begin{split}L(\theta)&=\prod_{j=1}^{n} p_j(\theta)^{x_j}\left(1-p_j(\theta)\right)^{1-x_j} =\prod_{j=1}^{n}\left(\dfrac{1}{ 1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\right)^{x_j} \left(\dfrac{\exp \left[ -Da\left(\theta -b_{j}\right) \right]}{ 1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\right)^{1-x_j}\\ &= \prod_{j=1}^{n} \dfrac{\left(\exp \left[ -Da\left(\theta -b_{j}\right) \right]\right)^{1-x_j}}{ 1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]} \end{split}\).

対数尤度関数

\(\begin{split} l(\theta)&=\sum_{j=1}^{n}\left[-Da(\theta -b_{j})(1-x_j)\right] -\sum_{j=1}^{n}\log\left(1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]\right) \end{split}\)

　　　\(\displaystyle\therefore \sum_{j=1}^{n}x_j =\sum_{j=1}^{n}\frac{1}{1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\).
項目反応パタン\((x_1,x_2,\cdots,x_n)\) に対し、この方程式を満たす\(\theta\) が能力パラメータということになる。
この方程式は代数的に解くことができないので、逐次近似法により近似解を求める。

1.1 ニュートン法.

方程式 \(f(x)=0\)　の近似解は.

　　　\({\begin{aligned} x_{n+1} &=x_n - \frac{f(x_n)}{f'(x_n)} \end{aligned} }\)　によって逐次近似で求まる。

対数尤度関数 \(l'(\theta)=0\) に対して適用すると、
　　　\(\theta_{n+1} =\theta_n - \dfrac{l'(\theta_n)}{l''(\theta_n)}\)　（式はモデルによって異なる　p159参照）。

最尤法は尤度関数・対数尤度関数が最大となる値を推定値とする方法である。最尤推定値は、尤度関数の導関数が０となる値であり、ニュートン法で求めることが可能になる。

1.2 Rで最尤推定(2PL).

表5.1（p102参照）.

a <-c(1,1.5,0.7,1.2)
b <-c(1.5,0.5,0,-1)
u <-c(1,0,1,1)
t0 <-log(sum(u)/(length(u)-sum(u)))#能力の初期値
t0

## [1] 1.098612

能力の初期値に対する各項目の正答確率.

項目特性関数(ICF).

icc2PL <- function(a, b, theta){
  1/(1+exp(-1.7*a*(theta-b)))
}
icc2PL(a,b,t0)

## [1] 0.3357350 0.8214880 0.7870691 0.9863621

対数尤度関数.

LL <- function(u,a,b,theta){ sum(u*log(icc2PL(a,b,theta)) + (1-u)*log(1-icc2PL(a,b,theta)))
}
LL(u,a,b,t0)

## [1] -3.067703

一次導関数.

dLL <- function(u,a,b,theta){ 1.7*sum(a*(u-icc2PL(a,b,theta)))
} 
dLL(u,a,b,t0)

## [1] -0.6843347

二次導関数.

ddLL <- function(a,b,theta){ -1.7^2*sum(a^2*icc2PL(a,b,theta) * (1-icc2PL(a,b,theta))) #
}
ddLL(a,b,t0)

## [1] -1.891389

推定値の更新(１回目).

t1 <- t0-dLL(u,a,b,t0)/ddLL(a,b,t0)
t1

## [1] 0.7367963

収束までのアルゴリズム.

e <- 0.0001 # 収束基準
t0 <-log(sum(u)/(length(u)-sum(u)))#能力の初期値
#LL(u,a,b,t0)##対数尤度

conv <- 0
while(conv==0){
  t1 <- t0-dLL(u,a,b,t0)/ddLL(a,b,t0)
  if(abs(t1-t0)<e) conv <-1
  t0 <-t1
}
t1

## [1] 0.776438

LL(u,a,b,t1)

## [1] -2.953549

最適解と対数尤度が出力される。

1.2.1 フィッシャー情報量.

スコア関数.

　\(U(\theta;x) = \frac{\partial}{\partial \theta} \log L(\theta; x) = \frac{1}{L(\theta;x)} \frac{\partial}{\partial \theta} L(\theta;x)\).

フィッシャー情報量. \(\ I_{X}(\theta) = Var[U(\theta ;x)] = E\left[\left( \frac{\partial}{\partial \theta} \log L(\theta; x) \right)^{2}\right] = - E\left[\frac{\partial ^{2}}{\partial \theta ^{2}} \log L(\theta; x)\right]\).
対数尤度関数の２次導関数の期待値であり、曲率や最大値に関係する。フィッシャー情報量の平方根の逆数は最尤推定値の標準誤差であり、これを用いて区間推定を行うことができる。

2 ベイズ推定法.

2.1 ベイズの発見.

ベイズの発見を紐解いてみましょう。（想像です）.
ベイズは牧師だったそうです。ベイズはLost Coinの寓話が好きで、いつもお気に入りのコインを投げて遊んでいました。
でも、ベイズお気に入りのコインは表が出やすいような気がします。
　　４回投げてみると、表が３回出ました。
弟のコインを借りて４回投げてみると、表が２回出ました。
弟のコインは、表の出る確率 \(p\) は1/2のような気がします。

2.2 尤度.

高校で学ぶ二項分布では、表の出る確率 \(p\) としたとき、

表が２回出る確率は \({}_{4} C_{2}\:p^2\:(1-p)^2\).
この式を \(p\) の関数と見ることができるのではないか、というのがベイズの発見です。

これは４次関数です。
微積分の好きなベイズは微分してグラフを描くことにしました。

2.2.1 弟のコイン.

curve(choose(4,2)*x^2*(1-x)^2,from=0 ,to=1)

\(p=0.5\) のとき表が２回出る確率は最大なので、 \(p=0.5\)が最もありそうな値として推定されることになります。
このような推定法を最尤法といいます。
ベイズの弟のコインで表の出る確率は、0.5 と判断して良さようです。

2.2.2 ベイズのコイン.

４回投げて３回表が出る確率は \({}_{4} C_{3}\:p^3\:(1-p)\).

curve(choose(4,3)*x^3*(1-x)^1,from=0 ,to=1)

\(L(p)={}_{4} C_{3}\:p^3\:(1-p)\) として、
　　　\(L'(p)={}_{4} C_{3}\:(3p^2-4p^3)=0\) から、 \(p=3/4=0.75\). 最尤推定値は微分を利用して求めることができます。

上のグラフは、\(X=p\) を確率変数と考えたときの確率密度関数と同等であり、尤度関数は\(p\) の確率分布を表していると考えることができます。
つまり、尤度関数を \(p\) の確率分布と考えるのが尤度の考え方です。

2.3 事前分布.

４回のうち３回表が出たからといって、\(p=3/4=0.75\)と判断するのでは、相対度数を考えることと変りありません。
表と裏は同様に確からしいのですから、\(p=0.5\)というのが事前の予想としてあります。例えば、\(p\sim N\left( 0.5 ,0.1 ^{2}\right)\).

curve(dnorm(x,0.5,0.1),0,1)

これを事前の予測と考えます。pの確率密度関数で、pの事前分布といいます。
予測を確率分布で考えるところがポイントです。

2.4 事後分布.

情報が入ると確率判断が変化します。

乗法定理　\(P(H\cap D)=P(H)\cdot P(D|H)=P(D)\cdot P(H|D)\)から.
　　　\(P(H|D)=\dfrac{P(H)\cdot P(D|H)}{P(D)}\)　（分母は正規化定数で確率を１にする）

データD取得のもとでの条件確率の最大を考えるには、
　　　\(P(H|D) \propto P(H) \cdot P(D|H)\)　の最大を考えればよいことになります。

ベイズコインの事前分布の例では、

\(P(H)\)は表が出る確率がpである確率を意味し、
　　　\(P(H) = \dfrac{1}{\sqrt{2 \pi}\cdot0.1} \exp \left[ -\dfrac{(p-0.5)^2}{2 (0.1)^2} \right]\)
\(P(D|H)\)は表が出る確率がpであるときDが起きる確率ということになり、
　　　\(P(D|H)={}_{4} C_{3}\:p^3\:(1-p)\)

従って、ベイズの定理は、定数倍を無視すると.

　　　\(P(H|D) \propto \frac{1}{\sqrt{2 \pi}\cdot0.1}\exp \left[ -\frac{(p-0.5)^2}{2 (0.1)^2} \right] \cdot p^3(1-p)\)　　

curve(dnorm(x, mean=0.5,sd=0.1)*x^3*(1-x)^1,0,1,ylab="P(D|H)",xlab="p")

事後分布が最大となるpは、尤度関数をかけた影響で少し右へよっていることが分かります。これがベイズ推定による推定値で、データD取得により、pの推定値が更新されたことになります。
導関数が０となるpを求めればよいのですが、逐次近似法で近似計算するのがよいでしょう。

数理統計では、事前分布として共役分布を用いると計算が簡単になることが知られています。

このように、得られた情報を推定値に反映させる方法がベイズ推定です。

2.5 ［ベイズ推定のまとめ］.

予測は確率分布で表される。(確率密度関数).
情報がない段階での予測を 事前分布 という.
情報が入ると予測は変化する.
（勝負に勝てば次も勝つのではないかと考える）.

　これを予測のベイズ更新という。更新された予測を 事後分布 という。

2.6 事後確率最大(MAP).

maxuimum a posteriori.
事後分布は取得したデータDによって更新されたpの確率分布であり、MAPでは事後確率を最大とするpが推定値となります。
pは極大点であることから、導関数を０とする解を問題に帰着されます。推定値はニュートン法を利用した逐次近似法で求められます。

3 まとめ.

前書きにある通り、この中級編は骨太で長く、特に項目パラメータの推定はp199にあるように難解です。EMアルゴリズムの理解も必要となるので、中断して次の章に移ります。

4 参考.

アルゴリズム

u.all <- matrix(c(0,0,0,1,
                  0,1,0,1,
                  1,1,1,0,
                  1,0,1,0,
                  1,0,1,1),ncol=4,byrow=TRUE)


ni <- nrow(u.all) # 受験者数 (I)
e <- 0.0001 # 収束基準
theta.mle <- numeric(ni) # 能力パラメタ推定値格納用

for(i in 1:ni){ # 受験者ごとにループ
  u <- u.all[i,] # i番目の受験者の反応ベクトル t0 <- log(sum(u)/(length(u)-sum(u))) # 初期値 conv <- 0 # 収束判定の結果(収束したら 1 とする) while(conv==0){ # 収束していない間は繰り返す
  t1 <- t0-dLL(u,a,b,t0)/ddLL(a,b,t0) # 更新 if(abs(t1-t0)<e) conv <- 1 # 収束を判定
  t0 <- t1 # 更新した解を保存
}
theta.mle[i] <- t1 # 最適解を保存 }
theta.mle # 推定値

## [1] 0.000000 0.000000 0.000000 0.000000 1.034176

p164.

第８章 パラメタ推定法