名義反応モデル

# データを生成
import numpy as np
import pandas as pd

def simulate_nrm(N=1000, J=10, K=4, seed=42):
    rng = np.random.default_rng(seed)
    theta = rng.normal(0, 1, size=N)

    # カテゴリ0を参照カテゴリ (a=0, b=0) とする
    a_free = np.sort(rng.normal(0, 1, size=(J, K - 1)), axis=1)
    b_free = rng.normal(0, 1, size=(J, K - 1))

    a_full = np.zeros((J, K))
    a_full[:, 1:] = a_free
    b_full = np.zeros((J, K))
    b_full[:, 1:] = b_free

    U = np.zeros((N, J), dtype=int)
    for i in range(N):
        for j in range(J):
            logits = a_full[j] * theta[i] + b_full[j]
            logits -= logits.max()
            P = np.exp(logits)
            P /= P.sum()
            U[i, j] = rng.choice(K, p=P)
    return U, theta, a_free, b_free

num_users = 1000
num_items = 20
U, true_theta, true_a, true_b = simulate_nrm(N=num_users, J=num_items)
df = pd.DataFrame(
    U,
    index=[f"user_{i+1}" for i in range(num_users)],
    columns=[f"question_{j+1}" for j in range(num_items)],
)
df.head()

# 縦持ちへ変換
df_long = pd.melt(
    df.reset_index(),
    id_vars="index",
    var_name="item",
    value_name="response",
).rename(columns={"index": "user"}).astype({"user": "category", "item": "category"})

df_long.head()

import pymc as pm
import pytensor.tensor as pt

user_idx = df_long["user"].cat.codes.to_numpy()
users = df_long["user"].cat.categories.to_numpy()
item_idx = df_long["item"].cat.codes.to_numpy()
items = df_long["item"].cat.categories.to_numpy()
responses = df_long["response"].to_numpy().astype("int64")

K = int(responses.max() + 1)
n_free = K - 1
n_items = len(items)

coords = {
    "user": users,
    "item": items,
    "category_free": np.arange(n_free),
    "obs_id": np.arange(len(df_long)),
}

with pm.Model(coords=coords) as model:
    user_idx_ = pm.Data("user_idx", user_idx, dims="obs_id")
    item_idx_ = pm.Data("item_idx", item_idx, dims="obs_id")

    theta = pm.Normal("theta", 0.0, 1.0, dims="user")
    a_free = pm.Normal("a_free", 0.0, 1.0, dims=("item", "category_free"))
    b_free = pm.Normal("b_free", 0.0, 2.0, dims=("item", "category_free"))

    # 参照カテゴリ（0列目）にゼロを付与
    a_full = pt.concatenate([pt.zeros((n_items, 1)), a_free], axis=1)
    b_full = pt.concatenate([pt.zeros((n_items, 1)), b_free], axis=1)

    # logits[n, k] = a_{j,k} * theta_i + b_{j,k}
    logits = a_full[item_idx_] * theta[user_idx_][:, None] + b_full[item_idx_]
    p = pm.math.softmax(logits, axis=1)

    pm.Categorical("obs", p=p, observed=responses, dims="obs_id")

%%time
with model:
    idata = pm.sample(random_seed=0, draws=1000)

post_mean = idata.posterior.mean(dim=["chain", "draw"])

import matplotlib.pyplot as plt

fig, axes = plt.subplots(figsize=[12, 4], ncols=3)

ax = axes[0]
ax.scatter(true_theta, post_mean["theta"])
ax.plot(true_theta, true_theta, color="gray")
_ = ax.set(xlabel="true_theta", ylabel="theta_hat")

ax = axes[1]
ax.scatter(true_a.flatten(), post_mean["a_free"].to_numpy().flatten())
ax.plot(true_a.flatten(), true_a.flatten(), color="gray")
_ = ax.set(xlabel="true_a", ylabel="a_hat")

ax = axes[2]
ax.scatter(true_b.flatten(), post_mean["b_free"].to_numpy().flatten())
ax.plot(true_b.flatten(), true_b.flatten(), color="gray")
_ = ax.set(xlabel="true_b", ylabel="b_hat")