Go to Qiita Advent Calendar 2024 Top
LoginSignup
2
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@takuyakubo(久保 卓也)

試験の数理 その4（問題paramter推定の実装）

Last updated at Posted at 2020-09-17

試験の数理 その3（3PL modelの最適化） の続きです。

前回は、「3PL modelの最適化の数理について」でした。
今回は、「3PL modelの問題paramter推定の実装について」です。

用いた環境は、

  • python 3.8
  • numpy 1.19.2
  • scipy 1.5.2

です。

データ生成

試験の数理 その１（問題設定とデータの生成）で行ったように実験のために各種パラメータや項目反応行列を生成します。

import numpy as np
from functools import partial

# 数値計算で発散させないように小さなεをおいておく
epsilon =  0.0001
# 3 parameter logistic model の定義
def L3P(a, b, c, x):
    return c + (1 - epsilon - c) / (1 + np.exp(-  a * (x - b)))

# 2 parameter logistic model の定義。処理の統一のためにcも引数に取ることとする。
def L2P(a, b, c, x):
    return (1 - epsilon) / (1 + np.exp(-  a * (x - b)))

# model parameterの定義
# aは正の実数, bは実数, cは0より大きく1未満であれば良い

a_min = 0.7
a_max = 4

b_min = -2
b_max = 2

c_min = 0
c_max = .4

# 何問、何人にするか、下なら10問4000人
num_items = 10
num_users = 4_000

# 問題parameterの生成
item_params = np.array(
    [np.random.uniform(a_min, a_max, num_items),
     np.random.uniform(b_min, b_max, num_items),
     np.random.uniform(c_min, c_max, num_items)]
).T

# 受験者parameterの生成
user_params = np.random.normal(size=num_users)

# 項目反応行列の作成、 要素は1(正答)か0(誤答)
# i行j列は問iに受験者jがどう反応したか
ir_matrix_ij = np.vectorize(int)(
    np.array(
        [partial(L3P, *ip)(user_params) > np.random.uniform(0, 1, num_users) for ip in item_params]
    )
)

今までと同様に、問題を表す添字として$i$、受験者を表す添字として$j$を用います。
受験者を1000人以上に大きくとります。これは、経験的に1PLの推定には100人以上、2PLの推定には300人以上、3PLの推定には1000人以上あるとある程度推定が安定することが知られているからです。気になる方はここを変更して実験してみてください。

ここでは、問題parameterとして

a / 識別力 b / 困難度 c / 当て推量
問 1 3.34998814 0.96567682 0.33289520
問 2 1.78741502 1.09887666 0.22340858
問 3 1.33657604 -0.97455532 0.21594273
問 4 1.05624284 0.84572140 0.11501424
問 5 1.21345944 1.24370213 0.32661421
問 6 3.22726757 -0.95479962 0.33023057
問 7 1.73090248 1.46742090 0.21991115
問 8 2.16403443 1.66529355 0.10403063
問 9 2.35283349 1.78746377 0.22301869
問 10 1.77976105 -0.06035497 0.29241184

を得ました。これを推定するのが目的となります。

データの前処理

推定を行う前に推定の精度を安定させるために、

  • 正答率が高すぎる問題
  • 正答率が低すぎる問題
  • 素点との相関が低い問題

を推定対象から除去します。具体的には

# 推定しづらい問題を除去する。
# 正答率が高すぎたり(> 0.9)、低すぎたり(0.1 <)、素点との相関があまりにもない(< 0.3)問題は外す
filter_a = ir_matrix_ij.sum(axis=1) / num_users < 0.9
filter_b = ir_matrix_ij.sum(axis=1) / num_users > 0.1
filter_c = np.corrcoef(np.vstack((ir_matrix_ij, ir_matrix_ij.sum(axis=0))))[num_items][:-1] >= 0.3
filter_total = filter_a & filter_b & filter_c

# 項目反応行列を再定義する
irm_ij = ir_matrix_ij[filter_total]
num_items, num_users = irm_ij.shape

推定用の関数と定数の準備

前回見たように問題$i$毎に

\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}

を計算するのがE stepであり、

\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}

を0にするのが、M stepになります。
ここで、周辺化や上の積分は区分求積的に解くこととして、$\theta$の代表点として$X_k$や$g_k = g(X_k)$を用意します。

そのため、次を用意します。

# 受験者パラメータの取りうる範囲を定義する。
X_k = np.linspace(-4, 4, 41)
# 受験者パラメータの分布を定義する。ここでは、scipyの正規分布を使う。
from scipy.stats import norm
g_k = norm.pdf(X_k)
# E stepの関数
def get_exp_params(irm_ij, g_k, P3_ik):
    Lg_jk = np.exp(irm_ij.T.dot(np.log(P3_ik)) + (1 - irm_ij).T.dot(np.log(1 - P3_ik)))* g_k
    n_Lg_jk = Lg_jk / Lg_jk.sum(axis=1)[:, np.newaxis]
    f_k = n_Lg_jk.sum(axis=0)
    r_ik = irm_ij.dot(n_Lg_jk)
    return f_k, r_ik
# M step用のスコア関数
def score_(param, f_k, r_k, X_k):
    a, b, c = param
    P3_k = partial(L3P, a, b, c)(X_k)
    P2_k = partial(L2P, a, b, c)(X_k)
    R_k = r_k / P3_k - f_k
    v = [
        ((X_k - b) * R_k * P2_k).sum(),
        - a * (R_k * P2_k).sum(),
        R_k.sum() / (1 - c)
    ]
    return np.linalg.norm(v)

推定の実行

では、推定を実行します。推定は初期parameterに依存し、安定性もそこまで良くないことが予想されるので、初期parameterをランダムにいくつか用意して、安定したものの中の中央値を推定結果として採用します。M stepの最適化はscipyのminimizeを用いることとします。

from scipy.optimize import minimize
# minimize用の制約条件を定義する。
cons_ = {
    'type': 'ineq',
    'fun': lambda x:[
        x[0] - a_min,
        a_max - x[0],
        x[1] - b_min,
        b_max - x[1],
        x[2] - c_min,
        c_max - x[2],
    ]
}
# 初期parameter生成用のparameterを用意する。
a_min, a_max = 0.1, 8.0
b_min, b_max = -4.0, 4.0
c_min, c_max = epsilon, 0.6

# 推定実行用のoarameter
# EM algorithmの繰り返し終了条件
delta = 0.001
# EM algorithmの繰り返し最大回数
max_num_of_itr = 1000

# 数値安定のために何度か計算して、安定したものの中の中央値を採用する
p_data = []
for n_try in range(10):
    # 推定の初期値を定義する。
    item_params_ = np.array(
        [np.random.uniform(a_min, a_max, num_items),
         np.random.uniform(b_min, b_max, num_items),
         np.random.uniform(c_min, c_max, num_items)]
    ).T

    prev_item_params_ = item_params_
    for itr in range(max_num_of_itr):
        # E step : exp paramの計算 
        P3_ik = np.array([partial(L3P, *ip)(X_k) for ip in item_params_])
        f_k, r_ik = get_exp_params(irm_ij, g_k, P3_ik)
        ip_n = []
        # 各問題ごとに最適化問題をとく
        for item_id in range(num_items):
            target = partial(score_, f_k=f_k, r_k=r_ik[item_id], X_k=X_k)
            result = minimize(target, x0=item_params_[item_id], constraints=cons_, method="slsqp")         
            ip_n.append(list(result.x))

        item_params_ = np.array(ip_n)
        # 前回との平均差分が一定値を下回ったら計算終了
        mean_diff = abs(prev_item_params_ - item_params_).sum() / item_params_.size
        if mean_diff < delta:
            break
        prev_item_params_ = item_params_

    p_data.append(item_params_)
    
p_data_ = np.array(p_data)
result_ = []
for idx in range(p_data_.shape[1]):
    t_ = np.array(p_data)[:, idx, :]
    # 計算結果で極端なものを排除
    filter_1 = t_[:, 1] < b_max - epsilon 
    filter_2 = t_[:, 1] > b_min + epsilon
    # 残った中のmedianを計算結果とする。
    result_.append(np.median(t_[filter_1 & filter_2], axis=0))
    
result = np.array(result_)

実験結果

今回は計算の結果次を得ました。
計算結果(目的)

a / 識別力 b / 困難度 c / 当て推量
問 1 3.49633348(3.34998814) 1.12766137(0.96567682) 0.35744497(0.33289520)
問 2 2.06354365(1.78741502) 1.03621881(1.09887666) 0.20507606(0.22340858)
問 3 1.64406087(1.33657604) -0.39145998(-0.97455532) 0.48094315(0.21594273)
問 4 1.47999466(1.05624284) 0.95923840(0.84572140) 0.18384673(0.11501424)
問 5 1.44474336(1.21345944) 1.12406269(1.24370213) 0.31475672(0.32661421)
問 6 3.91285332(3.22726757) -1.09218709(-0.95479962) 0.18379076(0.33023057)
問 7 1.44498535(1.73090248) 1.50705016(1.46742090) 0.20601461(0.21991115)
問 8 2.37497907(2.16403443) 1.61937999(1.66529355) 0.10503096(0.10403063)
問 9 3.10840278(2.35283349) 1.69962392(1.78746377) 0.22051818(0.22301869)
問 10 1.79969976(1.77976105) 0.06053145(-0.06035497) 0.29944448(0.29241184)

問3のb(困難度)など若干おかしな数値になっているところもありますが、おおよそ良い精度で推定できています。

次回

次回は問題のparameterの結果を用いて受験者の能力parameter $\theta_j$の推定を行います。
試験の数理 その5（受験者paramter推定）

参考文献

2
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?