一般化部分採点モデル

# ESC（単調ケース）
E_mono = np.sum((np.arange(1, K+1)[None, :]) * P_mono, axis=1)

# ESC（非単調ケース）
E_nonmono = np.sum((np.arange(1, K+1)[None, :]) * P_nonmono, axis=1)

plt.figure(figsize=(5, 3))
plt.plot(theta, E_mono, label="monotone b")
plt.plot(theta, E_nonmono, linestyle="--", label="non-monotone b")

plt.xlabel(r"$\theta$")
plt.ylabel(r"$\mathbb{E}[Y\mid\theta]$")
plt.title("Expected Score Curve (ESC)")
plt.legend()
plt.tight_layout()
plt.show()

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

def simulate_gpcm(N=1000, J=10, K=4, seed=42):
    rng = np.random.default_rng(seed)
    theta = rng.normal(0, 1, size=N)
    a = rng.lognormal(0, 0.3, size=J)
    b = rng.normal(0, 1, size=(J, K - 1))  # step difficulties（順序制約なし）

    U = np.zeros((N, J), dtype=int)
    for j in range(J):
        inc = a[j] * (theta[:, None] - b[j][None, :])  # (N, K-1)
        csum = np.cumsum(inc, axis=1)                    # (N, K-1)
        log_num = np.concatenate([np.zeros((N, 1)), csum], axis=1)  # (N, K)
        mx = log_num.max(axis=1, keepdims=True)
        exp_shift = np.exp(log_num - mx)
        P = exp_shift / exp_shift.sum(axis=1, keepdims=True)
        for i in range(N):
            U[i, j] = rng.choice(K, p=P[i])
    return U, theta, a, b

num_users = 1000
num_items = 20
U, true_theta, true_a, true_b = simulate_gpcm(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)],
)

# 縦持ちへ変換
user_categories = [f"user_{i+1}" for i in range(num_users)]
item_categories = [f"question_{j+1}" for j in range(num_items)]

df_long = pd.melt(
    df.reset_index(), id_vars="index", var_name="item", value_name="response",
).rename(columns={"index": "user"})
df_long["user"] = pd.Categorical(df_long["user"], categories=user_categories, ordered=True)
df_long["item"] = pd.Categorical(df_long["item"], categories=item_categories, ordered=True)

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_steps = K - 1

coords = {
    "user": users,
    "item": items,
    "step": np.arange(n_steps),
    "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")
    responses_ = pm.Data("responses", responses, dims="obs_id")

    theta = pm.Normal("theta", 0, 1, dims="user")
    a = pm.LogNormal("a", 0, 0.5, dims="item")
    b = pm.Normal("b", 0, 2, dims=("item", "step"))

    # GPCM: logit_k = sum_{c=0}^{k-1} a_j (theta_i - b_{jc}),  logit_0 = 0
    inc = a[item_idx_][:, None] * (theta[user_idx_][:, None] - b[item_idx_])
    csum = pt.cumsum(inc, axis=1)
    logits = pt.concatenate([pt.zeros_like(csum[:, :1]), csum], axis=1)

    p = pt.special.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)

import matplotlib.pyplot as plt

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

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, post_mean["a"])
ax.plot(true_a, true_a, color="gray")
_ = ax.set(xlabel="true_a", ylabel="a_hat")

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

plt.tight_layout()