書記が数学やるだけ#823 項目反応理論モデルの推定

import pymc3 as pm
import numpy as np

# データの準備
responses = np.array([[1, 1, 1, 1, 0],
                      [1, 1, 1, 0, 0],
                      [1, 1, 0, 0, 0],
                      [1, 0, 0, 0, 0]])

num_items = len(responses[0])
num_people = len(responses)

# モデルの構築
with pm.Model() as irt_model:
    # 能力パラメータを事前分布を指定して定義
    ability = pm.Normal('ability', mu=0, sd=1, shape=num_people)
    
    # 項目パラメータ（難易度および識別力）を事前分布を指定して定義
    item_difficulty = pm.Normal('item_difficulty', mu=0, sd=1, shape=num_items)
    item_discrimination = pm.HalfNormal('item_discrimination', sd=1, shape=num_items)
    
    # 項目特性曲線（Item Characteristic Curve, ICC）を定義
    theta = ability[:, np.newaxis] - item_difficulty
    logits = item_discrimination * theta
    
    # 応答確率のロジスティック関数
    p = pm.math.invlogit(logits)
    
    # データ生成分布を指定
    observed = pm.Bernoulli('observed', p=p, observed=responses)

# 推定
with irt_model:
    trace = pm.sample(1000, tune=1000)

# 推定結果の確認
pm.summary(trace)
pm.traceplot(trace)
pm.plot_posterior(trace)

#テスト情報関数を計算する関数
def compute_test_information(ability_values, item_difficulty, item_discrimination):
    theta = ability_values[:, np.newaxis] - item_difficulty
    logits = item_discrimination * theta
    p = 1 / (1 + np.exp(-logits))
    information = item_discrimination**2 * np.exp(-logits) / (1 + np.exp(-logits))**2
    test_information = np.sum(information, axis=1)
    return test_information

# モデルからサンプリングされたパラメータの事後分布からランダムサンプルを抽出
sampled_ability = trace['ability'][::10]
sampled_item_difficulty = trace['item_difficulty'][::10]
sampled_item_discrimination = trace['item_discrimination'][::10]

# テスト情報関数の計算
ability_values = np.linspace(-3, 3, 1000)
test_information_values = np.zeros((len(sampled_ability), len(ability_values)))

for i in range(len(sampled_ability)):
    test_information_values[i, :] = compute_test_information(
        ability_values, sampled_item_difficulty[i, :], sampled_item_discrimination[i, :]
    )

# テスト情報関数の可視化
plt.figure(figsize=(10, 6))
plt.plot(ability_values, np.mean(test_information_values, axis=0), label='Mean Test Information')
plt.fill_between(
    ability_values,
    np.percentile(test_information_values, 2.5, axis=0),
    np.percentile(test_information_values, 97.5, axis=0),
    alpha=0.3,
    label='95% CI'
)
plt.title('Test Information Function')
plt.xlabel('Ability')
plt.ylabel('Test Information')
plt.legend()
plt.show()

#項目情報関数を計算する関数
def compute_item_information_function(ability_values, item_difficulty, item_discrimination):
    theta = ability_values[:, np.newaxis] - item_difficulty
    logits = item_discrimination * theta
    information = item_discrimination**2 * np.exp(-logits) / (1 + np.exp(-logits))**2
    return np.sum(information, axis=1)

# 計算と可視化
ability_values = np.linspace(-3, 3, 1000)  # ポイントの数を増やす

plt.figure(figsize=(12, 8))

for i in range(num_items):
    plt.subplot(2, 3, i + 1)
    
    # 寸法を一致させるためにability_valuesを再計算
    computed_ability_values = np.linspace(-3, 3, 1000)
    
    plt.plot(computed_ability_values, compute_item_information_function(
        computed_ability_values, sampled_item_difficulty[:, i], sampled_item_discrimination[:, i]
    ), label=f'Item {i + 1}')
    
    plt.title(f'Item Information Function - Item {i + 1}')
    plt.xlabel('Ability')
    plt.ylabel('Information')
    plt.legend()

plt.tight_layout()
plt.show()