1"""Two-parameter logistic IRT (marginal maximum likelihood, EM with Gauss–Hermite quadrature).23Items are binary (0/1); missing entries contribute nothing to the likelihood (missing-at-random by design, no4imputation). Returns item discriminations a_j, difficulties b_j, EAP scores and posterior SDs for each unit.5Deterministic (no randomness). Written for transparency rather than speed (N ≈ 600 × J ≈ 30 is instant).6"""7from __future__ import annotations8import numpy as np9from scipy.special import expit101112def fit_2pl(X: np.ndarray, n_quad: int = 41, max_iter: int = 1000, tol: float = 1e-6, min_items: int = 3,13 prior_sd_log_a: float = 0.5, prior_sd_b: float = 2.0):14 """MAP variant: log-normal prior on a (mean 0, sd prior_sd_log_a) and normal prior on b (sd prior_sd_b) to keep15 items that are nearly constant in the sample from running to the boundary (Heywood-type cases)."""16 X = np.asarray(X, dtype=float)17 N, J = X.shape18 obs = ~np.isnan(X)19 Xz = np.where(obs, X, 0.0)20 nodes, weights = np.polynomial.hermite_e.hermegauss(n_quad) # standard normal prior21 weights = weights / weights.sum()22 a = np.ones(J); b = np.zeros(J)23 ll_old = -np.inf24 for it in range(max_iter):25 # E-step: posterior over nodes for each unit26 P = expit(a[None, :] * (nodes[:, None] - b[None, :])) # Q x J27 logP, log1P = np.log(P + 1e-300), np.log(1 - P + 1e-300)28 # log-lik of unit i at node q = sum_j obs_ij [x_ij logP_qj + (1-x_ij) log(1-P_qj)]29 L = Xz @ logP.T + ((1 - Xz) * obs) @ log1P.T # N x Q30 L += np.log(weights)[None, :]31 m = L.max(axis=1, keepdims=True)32 post = np.exp(L - m); norm = post.sum(axis=1, keepdims=True); post /= norm33 ll = float((m.squeeze() + np.log(norm.squeeze())).sum())34 # M-step: Newton updates per item on expected counts35 r = post.T @ Xz # Q x J expected correct36 n = post.T @ obs.astype(float) # Q x J expected tries37 for j in range(J):38 for _ in range(10):39 p = expit(a[j] * (nodes - b[j]))40 g_a = np.sum((r[:, j] - n[:, j] * p) * (nodes - b[j])) - (np.log(a[j]) / prior_sd_log_a ** 2) / a[j]41 g_b = np.sum((r[:, j] - n[:, j] * p) * (-a[j])) - b[j] / prior_sd_b ** 242 w = n[:, j] * p * (1 - p)43 H_aa = -np.sum(w * (nodes - b[j]) ** 2) - (1 - np.log(a[j])) / (prior_sd_log_a ** 2 * a[j] ** 2)44 H_bb = -np.sum(w * a[j] ** 2) - 1 / prior_sd_b ** 245 H_ab = np.sum(w * a[j] * (nodes - b[j]))46 H = np.array([[H_aa, H_ab], [H_ab, H_bb]]) - 1e-6 * np.eye(2)47 step = np.linalg.solve(H, np.array([g_a, g_b]))48 a[j] -= step[0]; b[j] -= step[1]49 a[j] = float(np.clip(a[j], 0.05, 8.0)); b[j] = float(np.clip(b[j], -6, 6))50 if np.abs(step).max() < 1e-6:51 break52 if abs(ll - ll_old) < tol:53 break54 ll_old = ll55 # final EAP56 P = expit(a[None, :] * (nodes[:, None] - b[None, :]))57 L = Xz @ np.log(P + 1e-300).T + ((1 - Xz) * obs) @ np.log(1 - P + 1e-300).T + np.log(weights)[None, :]58 m = L.max(axis=1, keepdims=True); post = np.exp(L - m); post /= post.sum(axis=1, keepdims=True)59 eap = post @ nodes60 psd = np.sqrt(np.maximum(post @ nodes ** 2 - eap ** 2, 0))61 n_items = obs.sum(axis=1)62 eap = np.where(n_items >= min_items, eap, np.nan); psd = np.where(n_items >= min_items, psd, np.nan)63 return {"a": a, "b": b, "eap": eap, "psd": psd, "loglik": ll, "iterations": it + 1, "n_items": n_items}64