Skip to main content Instalar con Codex o Claude Copia este prompt, pégalo en Codex, Claude u otro asistente, y deja que revise la página de la skill y la instale por ti.
Un comando directo omite el prompt de revisión. Revisa el origen antes de ejecutarlo.
npx skills add https://github.com/xjtulyc/awesome-rosetta-skills --skill learning-experimentEl comando permanece en una sola línea. Desplázate horizontalmente para revisarlo antes de copiarlo.
¿Prefieres una copia local? Descarga los archivos que SkillsMP tiene disponibles ahora.
Ocupaciones relacionadasSOC
Basado en la clasificación ocupacional SOC
| name | learning-experiment |
| description | Educational psychology experiment design and analysis covering IRT, growth modeling, A/B testing, and learning curve estimation for research. |
| tags | ["education-research","item-response-theory","learning-analytics","experimental-design","growth-modeling"] |
| version | 1.0.0 |
| authors | ["@xjtulyc"] |
| license | MIT |
| platforms | ["claude-code","codex","gemini-cli","cursor"] |
| dependencies | {"python":["numpy>=1.24","pandas>=2.0","scipy>=1.11","statsmodels>=0.14","scikit-learn>=1.3","matplotlib>=3.7"]} |
| last_updated | 2026-03-17 |
| status | stable |
Learning Experiment Design and Analysis
When to Use This Skill
Use this skill when you need to:
- Design randomized controlled trials for educational interventions
- Analyze item response data with IRT models (1PL/2PL/3PL)
- Estimate student growth models and learning trajectories
- Conduct A/B testing for educational technology features
- Fit learning curves (power law, exponential decay of errors)
- Measure knowledge retention (forgetting curves, spacing effects)
- Compute test reliability (Cronbach's alpha, inter-rater reliability)
Trigger keywords: item response theory, Rasch model, 2PL, 3PL, student growth, learning curve, forgetting curve, spaced repetition, educational A/B test, mastery learning, knowledge tracing, Bayesian knowledge tracing, DIF analysis, NAEP, PISA, CAT, adaptive testing, effect size, pre-post design.
Background & Key Concepts
Item Response Theory (IRT)
The 3-parameter logistic model (3PL) gives the probability that student $\theta$ answers item $i$ correctly:
$$P(X_i=1|\theta) = c_i + (1-c_i) \cdot \frac{e^{a_i(\theta - b_i)}}{1 + e^{a_i(\theta - b_i)}}$$
- $a_i$: discrimination parameter
- $b_i$: difficulty parameter (item location)
- $c_i$: guessing parameter (lower asymptote)
The 2PL sets $c_i=0$; the 1PL (Rasch) additionally fixes $a_i=1$.
Item Information Function
$$I_i(\theta) = a_i^2 \cdot \frac{(P_i(\theta) - c_i)^2}{(1-c_i)^2} \cdot \frac{(1-P_i(\theta))}{P_i(\theta)}$$
Total test information: $I(\theta) = \sum_i I_i(\theta)$. Standard error: $SE(\theta) = 1/\sqrt{I(\theta)}$.
Power Law of Learning (Newell & Rosenbloom 1981)
$$T_n = T_1 \cdot n^{-\alpha}$$
where $T_n$ is the time (or error rate) on trial $n$ and $\alpha$ is the learning rate. Log-linearized: $\ln T_n = \ln T_1 - \alpha \ln n$.
Bayesian Knowledge Tracing (BKT)
Hidden Markov model with states $K \in {0,1}$ (unlearned/learned):
$$P(K_{n+1}=1 | K_n=0) = p_T \quad \text{(transit)}$$
$$P(X=1 | K=1) = 1-p_S,\quad P(X=1|K=0) = p_G$$
Effect Size Measures
- Cohen's $d$: $d = (\bar{x}_1 - \bar{x}_2)/s_p$ (pooled SD)
- Glass's $\Delta$: $\Delta = (\bar{x}_1 - \bar{x}2)/s{\text{control}}$
- Hedges' $g$: small-sample corrected version of $d$
Environment Setup
pip install numpy>=1.24 pandas>=2.0 scipy>=1.11 statsmodels>=0.14 \
scikit-learn>=1.3 matplotlib>=3.7
numpy np
pandas pd
scipy.optimize opt
statsmodels.api sm
matplotlib.pyplot plt
()
import
as
import
as
import
as
import
as
import
as
print
"Learning experiment analysis environment ready"
Core Workflow
Step 1: IRT Model Estimation (2PL)
import numpy as np
import pandas as pd
import scipy.optimize as opt
import matplotlib.pyplot as plt
np.random.seed(42)
n_students = 500
n_items = 20
theta_true = np.random.normal(0, 1, n_students)
a_true = np.random.uniform(0.5, 2.5, n_items)
b_true = np.random.linspace(-2.5, 2.5, n_items)
def p2pl(theta, a, b):
"""2PL probability of correct response."""
return 1.0 / (1.0 + np.exp(-a * (theta[:, None] - b[None, :])))
P_matrix = p2pl(theta_true, a_true, b_true)
responses = (np.random.uniform(size=P_matrix.shape) < P_matrix).astype(int)
print(f"Response matrix shape: {responses.shape}")
print(f"Mean item difficulty (p-values): {responses.mean(axis=0).round(2)}")
print(f"Overall correct rate: {responses.mean():.3f}")
def neg_log_likelihood_2pl(params, responses):
"""Compute -2LL for 2PL model with fixed quadrature points for theta.
Args:
params: vector of length 2*n_items (a1..an, b1..bn)
responses: (n_students, n_items) binary matrix
Returns:
negative log-likelihood
"""
n_i = responses.shape[1]
a = np.abs(params[:n_i]) + 0.2
b = params[n_i:]
n_quad = 15
quad_pts, quad_wts = np.polynomial.hermite.hermgauss(n_quad)
theta_q = quad_pts * np.sqrt(2)
wts = quad_wts / np.sqrt(np.pi)
total_ll = 0.0
for s in range(len(responses)):
p_iq = 1.0 / (1.0 + np.exp(-a[None, :] * (theta_q[:, None] - b[None, :])))
r = responses[s]
log_lik_q = (r * np.log(np.clip(p_iq, 1e-9, 1 - 1e-9)) +
(1 - r) * np.log(np.clip(1 - p_iq, 1e-9, 1 - 1e-9))).sum(axis=1)
lik_q = np.exp(log_lik_q)
marginal_lik = (lik_q * wts).sum()
total_ll += np.log(max(marginal_lik, 1e-300))
return -total_ll
a_init = np.ones(n_items)
b_init = np.zeros(n_items)
params_init = np.concatenate([a_init, b_init])
print("\nFitting 2PL model (simplified MML)...")
subset_students = responses[:100]
result = opt.minimize(neg_log_likelihood_2pl, params_init,
args=(subset_students,),
method="L-BFGS-B",
options={"maxiter": 50, "ftol": 1e-4})
a_est = np.abs(result.x[:n_items]) + 0.2
b_est = result.x[n_items:]
print(f"Optimization converged: {result.success}")
print(f"Correlation(a_true, a_est): {np.corrcoef(a_true, a_est)[0,1]:.3f}")
print(f"Correlation(b_true, b_est): {np.corrcoef(b_true, b_est)[0,1]:.3f}")
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
theta_plot = np.linspace(-4, 4, 200)
colors = ["steelblue", "orange", "green", "red"]
for idx, item_idx in enumerate([0, 5, 10, 15]):
p_true = 1 / (1 + np.exp(-a_true[item_idx] * (theta_plot - b_true[item_idx])))
axes[0].plot(theta_plot, p_true, color=colors[idx],
label=f"Item {item_idx+1} (b={b_true[item_idx]:.1f})")
axes[0].set_xlabel("θ (Ability)"); axes[0].set_ylabel("P(correct)")
axes[0].set_title("Item Characteristic Curves")
axes[0].legend(fontsize=8)
axes[0].axhline(0.5, color="gray", ls="--", lw=0.8)
for idx, item_idx in enumerate([0, 5, 10, 15]):
a, b = a_true[item_idx], b_true[item_idx]
p = 1 / (1 + np.exp(-a * (theta_plot - b)))
info = a**2 * p * (1 - p)
axes[1].plot(theta_plot, info, color=colors[idx],
label=f"Item {item_idx+1} (a={a:.2f})")
axes[1].set_xlabel("θ"); axes[1].set_ylabel("Information")
axes[1].set_title("Item Information Functions")
axes[1].legend(fontsize=8)
P_all = 1 / (1 + np.exp(-a_true[:, None] * (theta_plot[None, :] - b_true[:, None])))
test_info = (a_true[:, None]**2 * P_all * (1 - P_all)).sum(axis=0)
se_theta = 1 / np.sqrt(np.maximum(test_info, 0.01))
axes[2].plot(theta_plot, test_info, "steelblue", lw=2, label="Test Information")
ax2r = axes[2].twinx()
ax2r.plot(theta_plot, se_theta, "r--", lw=1.5, label="SE(θ)")
ax2r.set_ylabel("Standard Error", color="red")
axes[2].set_xlabel("θ"); axes[2].set_ylabel("Information")
axes[2].set_title("Test Information & SE")
axes[2].legend(loc="upper left")
plt.tight_layout()
plt.savefig("irt_analysis.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: irt_analysis.png")
Step 2: Learning Curves and Forgetting Models
import numpy as np
import pandas as pd
import scipy.optimize as opt
import matplotlib.pyplot as plt
np.random.seed(42)
n_students = 30
n_trials = 20
trials = np.arange(1, n_trials + 1, dtype=float)
T1_true = np.random.uniform(20, 40, n_students)
alpha_true = np.random.uniform(0.2, 0.6, n_students)
times = np.zeros((n_students, n_trials))
for s in range(n_students):
times[s] = T1_true[s] * trials**(-alpha_true[s]) * np.exp(
np.random.normal(0, 0.05, n_trials))
def power_law(n, T1, alpha):
"""Power law of practice."""
return T1 * n**(-alpha)
mean_times = times.mean(axis=0)
popt, pcov = opt.curve_fit(power_law, trials, mean_times,
p0=[30, 0.3], bounds=([1, 0.01], [200, 2.0]))
T1_fit, alpha_fit = popt
print(f"Power law fit: T1 = {T1_fit:.2f}s, α = {alpha_fit:.3f}")
print(f"95% CI for α: [{alpha_fit - 1.96*np.sqrt(pcov[1,1]):.3f}, "
f"{alpha_fit + 1.96*np.sqrt(pcov[1,1]):.3f}]")
def exponential_learning(n, T_inf, T1, k):
"""Exponential learning model."""
return T_inf + (T1 - T_inf) * np.exp(-k * (n - 1))
popt_exp, _ = opt.curve_fit(exponential_learning, trials, mean_times,
p0=[5, 35, 0.1], bounds=([0, 1, 0], [50, 200, 5]))
print(f"Exponential fit: T_inf = {popt_exp[0]:.2f}s, T1 = {popt_exp[1]:.2f}s, "
f"k = {popt_exp[2]:.4f}")
def aic(n_data, k_params, sse):
"""Compute AIC."""
sigma2 = sse / n_data
return n_data * np.log(sigma2) + 2 * k_params
sse_power = ((mean_times - power_law(trials, *popt))**2).sum()
sse_exp = ((mean_times - exponential_learning(trials, *popt_exp))**2).sum()
print(f"\nAIC (power law): {aic(n_trials, 2, sse_power):.2f}")
print(f"AIC (exponential): {aic(n_trials, 3, sse_exp):.2f}")
def ebbinghaus_retention(t, S):
"""Ebbinghaus forgetting curve: R(t) = e^{-t/S}."""
return np.exp(-t / S)
retention_days = np.array([1, 2, 7, 14, 30], dtype=float)
true_S = 12.0
retention_obs = ebbinghaus_retention(retention_days, true_S) + np.random.normal(0, 0.03, 5)
retention_obs = np.clip(retention_obs, 0, 1)
S_fit, _ = opt.curve_fit(ebbinghaus_retention, retention_days, retention_obs,
p0=[10], bounds=([1], [100]))
print(f"\nEbbinghaus S (stability) = {S_fit[0]:.1f} days")
print(f"Predicted 50% retention at: {S_fit[0] * np.log(2):.1f} days")
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for s in range(min(10, n_students)):
axes[0].plot(trials, times[s], color="lightblue", alpha=0.5, lw=0.8)
axes[0].plot(trials, mean_times, "k-", lw=2, label="Mean")
axes[0].plot(trials, power_law(trials, *popt), "r--", lw=2,
label=f"Power Law (α={alpha_fit:.2f})")
axes[0].set_xlabel("Trial"); axes[0].set_ylabel("Response Time (s)")
axes[0].set_title("Learning Curves")
axes[0].legend()
axes[1].plot(np.log(trials), np.log(mean_times), "ko", label="Data")
axes[1].plot(np.log(trials), np.log(power_law(trials, *popt)), "r-",
label=f"Fit (slope={-alpha_fit:.2f})")
axes[1].set_xlabel("ln(Trial)"); axes[1].set_ylabel("ln(Time)")
axes[1].set_title("Power Law Verification (Log-Log)")
axes[1].legend()
t_cont = np.linspace(0, 35, 300)
axes[2].scatter(retention_days, retention_obs, color="red", s=80, zorder=5,
label="Observed")
axes[2].plot(t_cont, ebbinghaus_retention(t_cont, S_fit[0]), "steelblue", lw=2,
label=f"Ebbinghaus (S={S_fit[0]:.1f}d)")
axes[2].set_xlabel("Days Since Learning"); axes[2].set_ylabel("Retention (fraction)")
axes[2].set_title("Forgetting Curve")
axes[2].set_ylim(0, 1.1)
axes[2].legend()
plt.tight_layout()
plt.savefig("learning_curves.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: learning_curves.png")
Step 3: A/B Testing Educational Interventions
import numpy as np
import pandas as pd
import scipy.stats as stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
np.random.seed(42)
n = 200
n_t = n // 2
pre_test = np.random.normal(55, 15, n)
pre_test = np.clip(pre_test, 0, 100)
treat = np.concatenate([np.ones(n_t), np.zeros(n_t)])
np.random.shuffle(treat)
te_individual = 8 + np.random.normal(0, 3, n)
post_test = (pre_test * 0.7
+ treat * te_individual
+ np.random.normal(0, 10, n))
post_test = np.clip(post_test, 0, 100)
df = pd.DataFrame({"student_id": range(n), "treat": treat,
"pre_test": pre_test, "post_test": post_test})
df["gain"] = df["post_test"] - df["pre_test"]
t_gains = df[df["treat"] == 1]["gain"]
c_gains = df[df["treat"] == 0]["gain"]
t_stat, p_val = stats.ttest_ind(t_gains, c_gains)
print("=== A/B Test: Treatment vs. Control ===")
print(f"Mean gain (treatment): {t_gains.mean():.2f}")
print(f"Mean gain (control): {c_gains.mean():.2f}")
print(f"Difference: {t_gains.mean() - c_gains.mean():.2f}")
print(f"t-statistic: {t_stat:.3f}, p-value: {p_val:.4f}")
pooled_sd = np.sqrt((t_gains.var() + c_gains.var()) / 2)
cohens_d = (t_gains.mean() - c_gains.mean()) / pooled_sd
j = 1 - 3 / (4 * (n - 2) - 1)
hedges_g = cohens_d * j
print(f"Cohen's d: {cohens_d:.3f}")
print(f"Hedges' g: {hedges_g:.3f}")
X_ancova = sm.add_constant(df[["treat", "pre_test"]])
ancova = sm.OLS(df["post_test"], X_ancova).fit(cov_type="HC3")
print("\n=== ANCOVA (post ~ treat + pre) ===")
print(ancova.summary().tables[1])
print(f"\nANCOVA treatment effect: {ancova.params['treat']:.3f} "
f"(SE={ancova.bse['treat']:.3f})")
from scipy.stats import norm
alpha_level = 0.05
observed_d = cohens_d
se_d = np.sqrt(2/n + observed_d**2 / (2*n))
power_obs = 1 - stats.norm.cdf(stats.norm.ppf(1 - alpha_level/2) - observed_d / se_d)
print(f"\nPost-hoc power: {power_obs:.3f}")
df["baseline_group"] = pd.cut(df["pre_test"], bins=[0, 50, 100],
labels=["Low", "High"])
subgroup = df.groupby(["baseline_group", "treat"])["gain"].agg(["mean", "std", "count"])
print("\n=== Subgroup Analysis (by baseline) ===")
print(subgroup.round(2))
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
colors = {1: "steelblue", 0: "orange"}
labels = {1: "Treatment (n=100)", 0: "Control (n=100)"}
for t_val in [0, 1]:
mask = df["treat"] == t_val
axes[0].scatter(df.loc[mask, "pre_test"], df.loc[mask, "post_test"],
color=colors[t_val], alpha=0.5, s=20, label=labels[t_val])
axes[0].plot([0, 100], [0, 100], "k--", lw=0.8, label="No change")
axes[0].set_xlabel("Pre-test Score"); axes[0].set_ylabel("Post-test Score")
axes[0].set_title("Pre-Post Scores by Treatment")
axes[0].legend(fontsize=8)
axes[1].hist(c_gains, bins=20, alpha=0.6, color="orange", label="Control", density=True)
axes[1].hist(t_gains, bins=20, alpha=0.6, color="steelblue", label="Treatment", density=True)
axes[1].axvline(c_gains.mean(), color="orange", ls="--", lw=2)
axes[1].axvline(t_gains.mean(), color="steelblue", ls="--", lw=2)
axes[1].set_xlabel("Score Gain"); axes[1].set_ylabel("Density")
axes[1].set_title(f"Distribution of Gains\n(d={cohens_d:.2f})")
axes[1].legend()
axes[2].scatter(ancova.fittedvalues, df["post_test"], alpha=0.4, s=15,
c=df["treat"].map({0: "orange", 1: "steelblue"}))
axes[2].plot([20, 100], [20, 100], "k--")
axes[2].set_xlabel("ANCOVA Fitted"); axes[2].set_ylabel("Observed")
axes[2].set_title(f"ANCOVA Fit (R²={ancova.rsquared:.3f})")
plt.tight_layout()
plt.savefig("educational_ab_test.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: educational_ab_test.png")
Advanced Usage
Bayesian Knowledge Tracing (BKT)
import numpy as np
from scipy.optimize import minimize
def bkt_forward(responses, p_L0, p_T, p_S, p_G):
"""Forward pass of BKT (single skill, single student).
Args:
responses: binary list of responses (1=correct, 0=incorrect)
p_L0: prior probability of knowing skill at trial 1
p_T: probability of transitioning to learned state
p_S: slip probability (knows but answers wrong)
p_G: guess probability (doesn't know but answers correctly)
Returns:
p_Ln: posterior P(learned) at each trial
log_likelihood
"""
p_Ln = np.zeros(len(responses) + 1)
p_Ln[0] = p_L0
log_lik = 0.0
for n, r in enumerate(responses):
p_L = p_Ln[n]
p_correct = p_L * (1 - p_S) + (1 - p_L) * p_G
log_lik += np.log(p_correct + 1e-10) if r == 1 else np.log(1 - p_correct + 1e-10)
if r == 1:
p_L_given_r = p_L * (1 - p_S) / (p_correct + 1e-10)
else:
p_L_given_r = p_L * p_S / (1 - p_correct + 1e-10)
p_Ln[n + 1] = p_L_given_r + (1 - p_L_given_r) * p_T
return p_Ln[1:], log_lik
np.random.seed(42)
true_params = {"p_L0": 0.2, "p_T": 0.15, "p_S": 0.08, "p_G": 0.20}
n_trials = 25
p_L_true = np.zeros(n_trials + 1)
p_L_true[0] = true_params["p_L0"]
for i in range(n_trials):
if np.random.random() > p_L_true[i]:
p_L_true[i+1] = p_L_true[i] + (1 - p_L_true[i]) * true_params["p_T"]
else:
p_L_true[i+1] = p_L_true[i]
learned = p_L_true[1:] > 0.5
sim_responses = []
for is_learned in learned:
if is_learned:
sim_responses.append(1 if np.random.random() > true_params["p_S"] else 0)
else:
sim_responses.append(1 if np.random.random() < true_params["p_G"] else 0)
def neg_ll(params):
p_L0, p_T, p_S, p_G = params
if not all(0.01 < p < 0.99 for p in params):
return 1e10
_, ll = bkt_forward(sim_responses, p_L0, p_T, p_S, p_G)
return -ll
result = minimize(neg_ll, x0=[0.3, 0.2, 0.1, 0.25], method="Nelder-Mead")
p_L0_f, p_T_f, p_S_f, p_G_f = result.x
print(f"BKT MLE: L0={p_L0_f:.3f}, T={p_T_f:.3f}, S={p_S_f:.3f}, G={p_G_f:.3f}")
p_mastery, _ = bkt_forward(sim_responses, p_L0_f, p_T_f, p_S_f, p_G_f)
print(f"Final mastery probability: {p_mastery[-1]:.3f}")
print(f"Trial when mastery > 0.95: "
f"{np.argmax(p_mastery > 0.95) + 1 if (p_mastery > 0.95).any() else 'Not reached'}")
Differential Item Functioning (DIF) Detection
import numpy as np
import pandas as pd
from scipy.stats import chi2_contingency
def mantel_haenszel_dif(responses, group, ability_strata=5):
"""Mantel-Haenszel DIF statistic for each item.
Args:
responses: (n_students, n_items) binary matrix
group: (n_students,) binary array — 0=reference, 1=focal
ability_strata: number of ability level strata
Returns:
DataFrame with item DIF statistics
"""
n_items = responses.shape[1]
total_score = responses.sum(axis=1)
strata = pd.qcut(total_score, ability_strata, labels=False, duplicates="drop")
results = []
for item in range(n_items):
A = B = C = D = 0
for s in range(ability_strata):
mask_s = strata == s
if mask_s.sum() < 4:
continue
ref_mask = mask_s & (group == 0)
foc_mask = mask_s & (group == 1)
n_s = mask_s.sum()
a = (responses[foc_mask, item]).sum()
b = (responses[ref_mask, item]).sum()
c = foc_mask.sum() - a
d = ref_mask.sum() - b
A += a * d / n_s
B += b * c / n_s
C += (a + d) * (a + c) * (b + c) * (b + d) / (n_s**2 * (n_s - 1) + 1e-10)
alpha_mh = A / (B + 1e-10)
delta_mh = -2.35 * np.log(alpha_mh + 1e-10)
results.append({"item": item, "MH_alpha": alpha_mh, "delta_MH": delta_mh})
return pd.DataFrame(results)
np.random.seed(42)
n_s = 400
n_i = 15
theta = np.random.normal(0, 1, n_s)
group = np.random.binomial(1, 0.4, n_s)
theta[group == 1] -= 0.5
a_dif = np.ones(n_i)
b_dif = np.linspace(-2, 2, n_i)
b_dif[5] += 0.8
P_dif = 1 / (1 + np.exp(-a_dif[None, :] * (theta[:, None] - b_dif[None, :])))
resp_dif = (np.random.uniform(size=P_dif.shape) < P_dif).astype(int)
dif_results = mantel_haenszel_dif(resp_dif, group)
print("=== DIF Analysis (Mantel-Haenszel) ===")
print(dif_results.round(3).to_string())
flagged = dif_results[dif_results["delta_MH"].abs() > 1.0]
print(f"\nFlagged items (|ΔMH| > 1.0): {flagged['item'].tolist()}")
Troubleshooting
| Problem | Cause | Fix |
|---|
| IRT estimation diverges | Extreme item p-values (< 0.05 or > 0.95) | Remove extreme items; increase sample |
| Power law fit R² < 0.80 | Plateau or ceiling effects | Check for ceiling; use exponential model |
| BKT mastery oscillates | Slip > 0.5 or guess > 0.5 | Constrain params: p_S < 0.4, p_G < 0.4 |
| Negative Cohen's d | Score direction | Verify higher score = better; check treatment coding |
| ANCOVA violation | Non-parallel slopes | Test treatment × pre interaction; use moderated regression |
| Negative DIF delta | Focal group advantaged | Expected; classify DIF direction (A favors reference, B favors focal) |
External Resources
- Lord, F. M. (1980). Applications of Item Response Theory. Educational Testing Service.
- Embretson, S. E., & Reise, S. P. (2000). Item Response Theory for Psychologists. Erlbaum.
- Anderson, J. R., et al. (1995). Cognitive tutors: Lessons learned. JLS, 4(2).
- mirt R package documentation
- EdPsych Measurement MOOC materials
Examples
Example 1: Test Reliability Analysis
import numpy as np
import pandas as pd
def cronbach_alpha(responses):
"""Compute Cronbach's alpha reliability coefficient.
Args:
responses: (n_students, n_items) matrix of item scores
Returns:
alpha coefficient
"""
n_items = responses.shape[1]
item_variances = responses.var(axis=0, ddof=1)
total_variance = responses.sum(axis=1).var(ddof=1)
alpha = (n_items / (n_items - 1)) * (1 - item_variances.sum() / total_variance)
return alpha
def item_statistics(responses):
"""Compute classical test theory item statistics."""
total = responses.sum(axis=1)
stats = []
for j in range(responses.shape[1]):
p_val = responses[:, j].mean()
corrected_total = total - responses[:, j]
r_it = np.corrcoef(responses[:, j], corrected_total)[0, 1]
stats.append({"item": j+1, "p_value": p_val, "r_item_total": r_it})
return pd.DataFrame(stats)
np.random.seed(42)
n_s, n_i = 300, 25
theta_r = np.random.normal(0, 1, n_s)
b_r = np.random.normal(0, 1, n_i)
P_r = 1 / (1 + np.exp(-(theta_r[:, None] - b_r[None, :])))
resp_r = (np.random.uniform(size=P_r.shape) < P_r).astype(int)
alpha = cronbach_alpha(resp_r)
print(f"Cronbach's α = {alpha:.3f}")
if alpha >= 0.90: print("Excellent reliability")
elif alpha >= 0.80: print("Good reliability")
elif alpha >= 0.70: print("Acceptable reliability")
else: print("Poor reliability — revise items")
item_stats = item_statistics(resp_r)
print("\nItem Statistics (first 5):")
print(item_stats.head().round(3))
poor_items = item_stats[item_stats["r_item_total"] < 0.20]
print(f"\nPoor discriminating items (r < 0.20): {poor_items['item'].tolist()}")
Example 2: Spaced Repetition Optimization
import numpy as np
import matplotlib.pyplot as plt
def sm2_schedule(initial_easiness=2.5, n_reviews=10, min_ease=1.3):
"""Simulate SM-2 spaced repetition algorithm schedule.
Args:
initial_easiness: starting E-factor (2.5 by default)
n_reviews: number of review sessions
min_ease: minimum E-factor
Returns:
DataFrame with intervals and retention estimates
"""
import pandas as pd
records = []
interval = 1
easiness = initial_easiness
stability = 1.0
for rev in range(1, n_reviews + 1):
q = 4
easiness = max(min_ease, easiness + 0.1 - (5 - q) * (0.08 + (5 - q) * 0.02))
if rev == 1:
next_interval = 1
elif rev == 2:
next_interval = 6
else:
next_interval = round(interval * easiness)
stability_at_review = stability * easiness
retention = np.exp(-interval / (stability_at_review * 10))
records.append({"review": rev, "interval_days": interval,
"easiness": easiness, "retention": retention})
stability = stability_at_review
interval = next_interval
return pd.DataFrame(records)
schedule = sm2_schedule()
print("SM-2 Spaced Repetition Schedule:")
print(schedule.to_string(index=False))
print(f"\nTotal study sessions to review 10 times: {schedule['interval_days'].sum()} days")
