試験の数理その2（項目反応理論の数理モデル）の続きです。

前回は、「項目反応理論で使われる数理モデルについて」でした。
今回は、「3PL modelの最適化の数理について」です。

今回は、数学的なお話のみです。ここで得られた結果を使って次回最適化を行います。
個人的に最も気合が入っているのはここですが、難しいと感じた場合は読み飛ばしても良いかもしれません。

推定方法について

まず、ここで扱う3PL modelについて復習します。前回扱ったとおり3PL modelは、以下で定義されます。

\Pr\{U_{ij} = 1|\theta, a, b, c\} = c_i + \frac{1 - c_i}{1 + \exp(-a_i(\theta_j - b_i))}

ここで、問題に関するparameterは$i$で、受験者に関するparameterは$j$で添字づけます。$U_{ij}$は問題$i$を受験者$j$が正答できるかどうかを示す2値確率変数で、1の時正答、0の時誤答です。各パラメータについては、

$a_i$は問題$i$の識別力で、とりうる範囲は正の実数
$b_i$は問題$i$の困難度で、とりうる範囲は実数全体
$c_i$は問題$i$の当て推量で、とりうる範囲は$[0, 1]$
$\theta_j$は受験者$j$の能力値で、とりうる範囲は実数全体

です。
いま、各受験者にとって、各問題は独立であると仮定し、

\begin{align}
P_{ij} &= \Pr\{U_{ij} = 1|\theta, a, b, c\}\\
Q_{ij} &= \Pr\{U_{ij} = 0|\theta, a, b, c\} = 1 - P_{ij}
\end{align}

とおくと、項目反応行列 $u = [u_{ij}]_{ij}$(ただし、$u_{ij} \in \{0, 1\}$)がえられる確率は、

\begin{align}
\Pr\{U = u|\theta, a, b, c\} &= \prod_{i, j} \Pr\{U_{ij} = u_{ij}|\theta, a, b, c\}\\
&= \prod_{i, j}P_{ij}^{u_{ij}}Q_{ij}^{1-u_{ij}}
\end{align}

となります。ここで、項目反応行列 $u$が与えられたときに、上の式を最大化するようなparameterを求めるのが、modelの最適化となるのですが、手法によって大きく３通りのやり方があります。

最尤推定法
ベイズ推定法
サンプリングによる方法(MCMCなど)

今回は、最尤推定法である周辺最尤推定法、中でも計算量もそこまで膨大にならない Bock & Aitkinの方法について述べます。周辺最尤推定法は受験者parameter $\theta$を周辺化して、まず先に問題parameter $a, b, c$を求める方法です。その中でも、Bock & Aitkinの方法は、受験者parameter $\theta$を隠れ変数とみてEMアルゴリズムを適用する方法となります。

記号の準備と簡単な計算

次の節で、Bock & Aitkinの方法についての具体的な計算をしますが、その前に計算に使う記号の準備とごく簡単な計算をここで準備しておきます。

導入する記号

$P_{ij}^*$: 2PL modelの式。つまり、 $P_{ij}^* = (1 + \exp(-a_i(\theta_j - b_i))^{-1}$
$Q_{ij}^* = 1 - P_{ij}^*$
$\partial^a_i = \partial^a = \frac{\partial}{\partial a_i}$
$\partial^b_i = \partial^b =\frac{\partial}{\partial b_i}$
$\partial^c_i = \partial^c =\frac{\partial}{\partial c_i}$
$\partial$: $\partial^a, \partial^b,\partial^c$のいずれか
$\partial^*$: $\partial^a, \partial^b$のいずれか

簡単な計算

$P_{ij} = c_i + (1 - c_i)P_{ij}^*$
$Q_{ij} = 1 - P_{ij} = (1 - c_i)(1 - P_{ij}^*) = (1 - c_i)Q_{ij}^* $
$\partial (\ln P_{ij}) = \partial P_{ij} / P_{ij}$
$\partial Q_{ij} = -\partial P_{ij}$
$\partial^* Q_{ij}^* = -\partial^* P_{ij}^*$
$\partial^* P_{ij} = (1 - c_i)\partial^* P_{ij}^*$
$\partial^c P_{ij} = 1 - P_{ij}^* = Q_{ij}^*$
$\partial^a P_{ij}^* = (\theta_j - b_i)P_{ij}^*Q_{ij}^*$
$\partial^b P_{ij}^* = -a_i P_{ij}^*Q_{ij}^*$

ここで、$\ln$は自然対数です。

Bock & Aitkinの方法

ここでは、問題parameterの推定方法としてBock & Aitkinの方法を述べます。Bock & Aitkinの方法は受験者parameter $\theta$を隠れ変数とみてEMアルゴリズムを適用する方法です。

EM アルゴリズム

EMアルゴリズムは観測量$x$, 非観測量$y$, モデルparameter $z$があったときに$p(x|y, z)$あるいは$\ln p(x|y, z)$を最大化させる$z$を探索するアルゴリズムになります。

具体的な手順は次のとおりです。

前提
$p(x, y | z)$が与えられているとする。

手順
1.　$z$を初期化する。次の手順も踏まえて$z^{old}$とする。
2.　(E step) $p(y | x, z^{old})$を計算する。
3.　(M step)　$z^{new} = \arg \max_z \sum_y p(y | x, z^{old}) \ln p(x, y | z) $　で更新する。
4.　収束を判定し、2, 3を繰り返す。

項目反応理論への適用

Bock & Aitkinの方法は受験者parameter $\theta$を隠れ変数とみてEMアルゴリズムを適用するとしましたが、今与えられているのは$\Pr\{U = u|\theta, a, b, c\} $だけです。ここで各受験者の能力$\theta_j$は同一の分布から独立に生成されるという仮定をいれると、その分布　$g(\theta_j)$をもちいて、

\begin{align}
\Pr\{U = u, \theta| a, b, c\} &= \prod_{j}\Pr\{U_j = u_j,|\theta_j, a, b, c\}g(\theta_j)\\
&= \prod_{i, j} P_{ij}^{u_{ij}}Q_{ij}^{1-u_{ij}}g(\theta_j)
\end{align}

と書けます。ところで、M stepについてここで最大化する量を書き下すと

\begin{align}
\mathcal{L}(a, b, c) &= \int \Pr\{\theta| U = u, a^{old}, b^{old}, c^{old}\} \ln\Pr\{U = u, \theta| a, b, c\}d\theta\\
&= \int \Pr\{\theta| U = u, a^{old}, b^{old}, c^{old}\} \sum_{i,j} u_{ij}\ln P_{ij} + (1 - u_{ij})\ln Q_{ij} + g(\theta_j) d\theta
\end{align}

となります。$\theta_j$の独立の仮定と、M stepでの最適化に関わる量$a, b, c$だけに関わるところだけに注目すると、最適化する量は

\begin{align}
\mathcal{\tilde{L}}(a, b, c) 
&= \sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} \sum_{i}(u_{ij}\ln P_{ij} + (1 - u_{ij})\ln Q_{ij} )d\theta_j
\end{align}

となります。ここで、$\mathcal{\tilde{L}}(a, b, c) $の最大にするparameterを探すわけですが、最大化されるということは、各parameterで偏微分 $\partial \mathcal{\tilde{L}}$ が0になることを意味しています。積分が微分と交換するとして、各微分を計算すると問題paramter $i$について

\begin{align}
\partial\mathcal{\tilde{L}}(a, b, c) 
&= \sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij}\partial\ln P_{ij} + (1 - u_{ij})\partial\ln Q_{ij} )d\theta_j\\
&= \sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij}\frac{\partial P_{ij}}{P_{ij}} + (1 - u_{ij})\frac{\partial Q_{ij}}{Q_{ij}} )d\theta_j\\
&= \sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} \frac{u_{ij} - P_{ij}}{P_{ij}Q_{ij}}\partial P_{ij}d\theta_j
\end{align}

よって各$a_i, b_i, c_i$について

\begin{align}
\partial^a\mathcal{\tilde{L}}(a, b, c) 
&= (1 - c_i)\sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij} - P_{ij})(\theta_j - b_i)\frac{P_{ij}^*Q_{ij}^*}{P_{ij}Q_{ij}}d\theta_j\\
&= \sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij} - P_{ij})(\theta_j - b_i)\frac{P_{ij}^*}{P_{ij}}d\theta_j\\
\partial^b\mathcal{\tilde{L}}(a, b, c) 
&= -a_i(1 - c_i)\sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij} - P_{ij})\frac{P_{ij}^*Q_{ij}^*}{P_{ij}Q_{ij}}d\theta_j\\
&= -a_i\sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij} - P_{ij})\frac{P_{ij}^*}{P_{ij}}d\theta_j\\
\partial^c\mathcal{\tilde{L}}(a, b, c) &= \sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\} \frac{u_{ij} - P_{ij}}{P_{ij}Q_{ij}} Q_{ij}^*d\theta_j \\
&= \frac{1}{(1 - c_i)}\sum_{j} \int \Pr\{\theta_j| U_j = u_j, a^{old}, b^{old}, c^{old}\}(u_{ij} - P_{ij}) \frac{1}{P_{ij}} d\theta_j
\end{align}

となります。ここで、$\theta_j$の独立同時性と$P_{ij}, P_{ij}^*$が$j$に対して、$\theta_j$のみで依存していることを考えると、
$\theta_j \rightarrow \theta, P_{ij}^{(*)} \rightarrow P_{i\theta}^{(*)}$として、

\begin{align}
\partial^a\mathcal{\tilde{L}}(a, b, c) 
&= \sum_{j} \int \Pr\{\theta| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij} - P_{i\theta})(\theta - b_i)\frac{P_{i\theta}^*}{P_{i\theta}}d\theta\\
\partial^b\mathcal{\tilde{L}}(a, b, c) 
&= -a_i\sum_{j} \int \Pr\{\theta| U_j = u_j, a^{old}, b^{old}, c^{old}\} (u_{ij} - P_{i\theta})\frac{P_{i\theta}^*}{P_{i\theta}}d\theta\\
\partial^c\mathcal{\tilde{L}}(a, b, c) 
&= \frac{1}{(1 - c_i)}\sum_{j} \int \Pr\{\theta| U_j = u_j, a^{old}, b^{old}, c^{old}\}(u_{ij} - P_{i\theta}) \frac{1}{P_{i\theta}} d\theta
\end{align}

ここで、

\begin{align}
r_{i\theta} &:= 
\sum_{j}u_{ij} \Pr\{\theta| U_j = u_j, a^{old}, b^{old}, c^{old}\} \\
f_\theta &= 
\sum_{j} \Pr\{\theta| U_j = u_j, a^{old}, b^{old}, c^{old}\} 
\end{align}

とおけば、

\begin{align}
\partial^a\mathcal{\tilde{L}}(a, b, c) 
&= \int  (r_{i\theta} - P_{i\theta}f_\theta)(\theta - b_i)\frac{P_{i\theta}^*}{P_{i\theta}}d\theta\\
\partial^b\mathcal{\tilde{L}}(a, b, c) 
&= -a_i \int (r_{i\theta} - P_{i\theta}f_\theta)\frac{P_{i\theta}^*}{P_{i\theta}}d\theta\\
\partial^c\mathcal{\tilde{L}}(a, b, c) 
&= \frac{1}{(1 - c_i)} \int (r_{i\theta} - P_{i\theta}f_\theta) \frac{1}{P_{i\theta}} d\theta
\end{align}

となり、さらに

\begin{align}
R_{i\theta} :=\frac{r_{i\theta}}{P_{i\theta}} - f_\theta
\end{align}

とおくと、

\begin{align}
\partial^a\mathcal{\tilde{L}}(a, b, c) 
&= \int R_{i\theta}(\theta - b_i)P_{i\theta}^*d\theta\\
\partial^b\mathcal{\tilde{L}}(a, b, c) 
&= -a_i \int R_{i\theta}P_{i\theta}^*d\theta\\
\partial^c\mathcal{\tilde{L}}(a, b, c) 
&= \frac{1}{(1 - c_i)} \int R_{i\theta} d\theta
\end{align}

となります。Bock & Aitkinの方法をまとめると、つぎのようになります。

$a, b, c$を適当に定める。それを$a^{old}, b^{old}, c^{old}$とする。
(E step +α) $R_{i\theta}$ をもとめる。
(M step) 各$i$について、$\partial\mathcal{\tilde{L}} = 0$となるような、$a^{new}, b^{new}, c^{new}$を求めて更新する。
4. Newton法や、Fissher point法など何らかの方法をもちいる。
収束を判定し、2, 3を繰り返す。

次回

今回の結果を踏まえて、python, numpyをつかって問題paramterの推定を実装します。
試験の数理その4（問題paramter推定の実装）

参考文献