Item response theory

Deep dive

How it works

In: binary answers x[i][j] (N persons, J items), or item parameters plus one person's answers · Out: (a_j, b_j, c_j) per item; θ̂ and SE per person · State: item table, fixed grid of K ability nodes · Cost: O(N·J·K) per EM iteration, O(J·K) to score a person

1. Model: P_j(θ) = c_j + (1 − c_j)·σ(a_j(θ − b_j)). Rasch (1PL): a_j = 1, c = 0. 2PL frees a_j; 3PL adds c_j, the floor at low ability (guessing). Local independence: given θ, answers are independent, so P(x | θ) = Π P_j^x (1−P_j)^(1−x).
2. Identify the scale: θ → s·θ + t, b → s·b + t, a → a/s leaves every probability unchanged. Fix it with a prior θ ~ N(0, 1) or by constraining the mean of b.
3. Fit items by marginal ML (Bock and Aitkin 1981): integrate θ out over the grid. E-step: posterior weights over the nodes per person. M-step: per item, a logistic regression of expected counts on node values, logit = aθ + c, concave so Newton converges (logistic regression, Newton). Repeat.
4. Score a person: EAP (posterior mean, posterior sd as SE), or MLE from Σ a_j(x_j − P_j) = 0 by Newton.
5. Information: 2PL I(θ) = Σ a_j² P_j Q_j; 3PL per item a²·(Q/P)·((P−c)/(1−c))². SE = 1/√I for the MLE; a N(0,1) prior adds about 1 to I.
6. Adaptive test: start at θ̂ = 0, ask the unused item with maximum I(θ̂), rescore, stop at the target SE.

Code

// 2PL: EM fit on a grid, EAP scoring, information, adaptive selection.
struct Rng(u64);
impl Rng {
    fn unit(&mut self) -> f64 { // splitmix64 -> (0,1)
        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        (((z ^ (z >> 31)) >> 11) as f64 + 0.5) / (1u64 << 53) as f64
    }
    fn normal(&mut self) -> f64 { (-2.0 * self.unit().ln()).sqrt() * (6.283185307179586 * self.unit()).cos() }
}
const K: usize = 41;
fn nodes() -> [f64; K] { std::array::from_fn(|k| -4.0 + 0.2 * k as f64) }
fn p(a: f64, b: f64, t: f64) -> f64 { 1.0 / (1.0 + (-a * (t - b)).exp()) }
/// Posterior over the grid under a N(0,1) prior (this prior fixes the ability scale).
fn post(it: impl Iterator<Item = (f64, f64, bool)>) -> [f64; K] {
    let th = nodes();
    let mut l: [f64; K] = std::array::from_fn(|k| -0.5 * th[k] * th[k]);
    for (a, b, x) in it {
        for k in 0..K { let q = p(a, b, th[k]); l[k] += if x { q.ln() } else { (1.0 - q).ln() }; }
    }
    let m = l.iter().cloned().fold(f64::MIN, f64::max);
    let s: f64 = l.iter().map(|v| (v - m).exp()).sum();
    std::array::from_fn(|k| (l[k] - m).exp() / s)
}
fn eap(w: &[f64; K]) -> (f64, f64) { // posterior mean and sd
    let th = nodes();
    let m: f64 = (0..K).map(|k| w[k] * th[k]).sum();
    (m, (0..K).map(|k| w[k] * (th[k] - m).powi(2)).sum::<f64>().sqrt())
}
fn info(a: f64, b: f64, t: f64) -> f64 { let q = p(a, b, t); a * a * q * (1.0 - q) } // Fisher information
/// Marginal ML by EM (Bock & Aitkin): E-step = posterior per person, M-step = weighted logistic fit per item.
fn fit(x: &[Vec<bool>], items: usize, iters: usize) -> Vec<(f64, f64)> {
    let th = nodes();
    let mut par = vec![(1.0, 0.0); items]; // (a, b)
    for _ in 0..iters {
        let (mut n, mut r) = (vec![[0.0; K]; items], vec![[0.0; K]; items]);
        for row in x {
            let w = post(row.iter().enumerate().map(|(j, &c)| (par[j].0, par[j].1, c)));
            for j in 0..items { for k in 0..K { n[j][k] += w[k]; if row[j] { r[j][k] += w[k]; } } }
        }
        for j in 0..items {
            let (mut a, mut c) = (par[j].0, -par[j].0 * par[j].1); // logit = a*t + c is concave: Newton is safe
            for _ in 0..4 {
                let (mut g, mut h) = ([0.0; 2], [0.0; 3]);
                for k in 0..K {
                    let q = 1.0 / (1.0 + (-(a * th[k] + c)).exp());
                    let (d, v) = (r[j][k] - n[j][k] * q, n[j][k] * q * (1.0 - q));
                    g[0] += d * th[k]; g[1] += d;
                    h[0] += v * th[k] * th[k]; h[1] += v * th[k]; h[2] += v;
                }
                let det = h[0] * h[2] - h[1] * h[1] + 1e-9;
                a = (a + (h[2] * g[0] - h[1] * g[1]) / det).clamp(0.1, 5.0);
                c += (h[0] * g[1] - h[1] * g[0]) / det;
            }
            par[j] = (a, -c / a);
        }
    }
    par
}
fn corr(u: &[f64], v: &[f64]) -> f64 {
    let n = u.len() as f64;
    let (mu, mv) = (u.iter().sum::<f64>() / n, v.iter().sum::<f64>() / n);
    let (mut s, mut su, mut sv) = (0.0, 0.0, 0.0);
    for (a, b) in u.iter().zip(v) { s += (a - mu) * (b - mv); su += (a - mu).powi(2); sv += (b - mv).powi(2); }
    s / (su * sv).sqrt()
}
fn main() {
    let mut rng = Rng(7);
    let (n, j) = (2000, 30);
    let truth: Vec<(f64, f64)> = (0..j).map(|_| ((0.3 * rng.normal()).exp(), rng.normal())).collect();
    let ability: Vec<f64> = (0..n).map(|_| rng.normal()).collect();
    let x: Vec<Vec<bool>> = ability.iter().map(|&t| truth.iter().map(|&(a, b)| rng.unit() < p(a, b, t)).collect()).collect();
    let est = fit(&x, j, 30);
    let (bt, be): (Vec<f64>, Vec<f64>) = (truth.iter().map(|t| t.1).collect(), est.iter().map(|t| t.1).collect());
    println!("corr difficulty {:.4}", corr(&bt, &be));
    // Posterior sd against the actual EAP error, for persons with |theta| < 1.
    let (mut sq, mut sd, mut cnt) = (0.0, 0.0, 0.0);
    for (row, &t) in x.iter().zip(&ability) {
        let (m, s) = eap(&post(row.iter().enumerate().map(|(i, &c)| (est[i].0, est[i].1, c))));
        if t.abs() < 1.0 { sq += (m - t).powi(2); sd += s * s; cnt += 1.0; }
    }
    println!("EAP rmse {:.3} vs mean posterior sd {:.3}", (sq / cnt).sqrt(), (sd / cnt).sqrt());
    let i0: f64 = truth.iter().map(|&(a, b)| info(a, b, 0.0)).sum();
    println!("test information at 0: {i0:.2}, 1/sqrt(I) = {:.3}", 1.0 / i0.sqrt());
    // CAT: 300-item bank, stop at posterior sd 0.3; max-information vs random order.
    let bank: Vec<(f64, f64)> = (0..300).map(|_| ((0.3 * rng.normal()).exp(), 1.2 * rng.normal())).collect();
    let (mut used, m) = ([0.0; 2], 200);
    for e in 0..m {
        for (mode, count) in used.iter_mut().enumerate() { // mode 0 = adaptive, 1 = random order
            let (mut left, mut got, mut th) = ((0..bank.len()).collect::<Vec<_>>(), vec![], 0.0);
            loop {
                let pick = if mode == 0 {
                    let i = |u: usize| info(bank[left[u]].0, bank[left[u]].1, th);
                    (0..left.len()).max_by(|&u, &v| i(u).total_cmp(&i(v))).unwrap()
                } else { (rng.unit() * left.len() as f64) as usize };
                let (a, b) = bank[left.swap_remove(pick)];
                got.push((a, b, rng.unit() < p(a, b, ability[e])));
                let (mean, sd) = eap(&post(got.iter().cloned()));
                th = mean;
                if sd < 0.3 || left.is_empty() { *count += got.len() as f64; break; }
            }
        }
    }
    println!("items to reach sd 0.3: adaptive {:.1}, random {:.1}", used[0] / m as f64, used[1] / m as f64);
    assert!(corr(&bt, &be) > 0.98);
    assert!(used[0] < used[1]);
}

Output, seed 7 (measured on this Mac (Apple M4 Pro)): difficulty correlation 0.9978; EAP rmse 0.343 for true |θ| < 1, mean posterior sd 0.356, 1/√I(0) = 0.359; adaptive 20.2 items, random 56.0. Seeds 1 to 10: correlation 0.993 to 0.998, adaptive 19.2 to 24.0 items, random 54.8 to 61.2; posterior sd 1 to 9 % above the rmse in all ten. Simulated abilities are normal, matching the prior.

Complexity and constants

• Fit O(N·J·K) per iteration, state 2·J·K floats; drivers are the exp and ln per response and node. Measured on this Mac (Apple M4 Pro), one thread, 30 iterations on 60,000 responses: 0.44 s at K = 41, 0.21 s at K = 21, 0.78 s at K = 81; correlation and rmse agreed to four decimals.

Efficiency

• Cost: 14.5 ms per EM iteration for 60,000 responses, 0.24 µs each (measured on this Mac (Apple M4 Pro)); compute-bound. Scoring one person from 20 answers: about 5 µs (one timing, same Mac).
• Levers:
  • Fewer grid nodes — K 41 → 21 halved fit time, same estimates (measured) — when ability is roughly normal.
  • Adaptive testing — 20.2 against 56.0 items for SE 0.3 (measured) — when the bank is large.
  • Rasch sum score — the likelihood depends on answers only through Σx (derivation), so scoring is a lookup of J+1 values — when all slopes are equal.
• Naive → standard → state of the art:
  • percent correct · O(N·J), confounds ability and difficulty
  • marginal ML by EM on a grid · O(N·J·K) per iteration
  • stochastic or variational fits, many dimensions · mirt (MHRM), py-irt

Implementation notes

• Posterior in log space (products of 100 probabilities underflow); the code clamps a to 0.1 to 5. With a = ln 10 / 400 per rating point, Rasch is the Elo expected score.

Pitfalls

• Joint ML (persons and items together) is biased at fixed test length: each person adds a parameter (Neyman and Scott 1948). Measured on this Mac (Apple M4 Pro), Rasch, 5,000 persons: difficulty spread 1.125 times too large at J = 10, 1.047 at J = 20, 1.018 at J = 40, near J/(J−1). Marginal ML avoids it.
• Perfect scores. All-correct (all-wrong) makes the likelihood rise (fall) without bound: MLE is +∞ (−∞). Use EAP, MAP or a censored score; EAP shrinks toward the prior mean, like the Bayesian average.
• Local independence breaks for items sharing a passage: information is overstated; group them as a testlet. One θ over mixed skills blends them.
• 3PL guessing: c shows only where P ≈ c, far below b; with few low-ability respondents it trades off against a and b. Fix it or use a prior.
• Differential item functioning: an item harder for one group at equal θ biases comparisons. Fit per group and test the difference (hypothesis testing).
• Adaptive selection favours b ≈ θ̂ and large a, overusing a few items: limit exposure, keep the bank large.

Variants and alternatives

In production

• Language exam: the Duolingo English Test is computer-adaptive, with item difficulty estimated from linguistic models (Settles et al. 2020).
• Health questionnaires: PROMIS adaptive tests pick bank items by earlier answers, 4 to 12 for the standard form (HealthMeasures; model not named).
• Model evaluation: IRT-built evaluation scales for NLP models (Lalor et al. 2016).

Implementations

No established Rust IRT crate found on crates.io.

Sources

• Georg Rasch, Probabilistic Models for Some Intelligence and Attainment Tests, 1960; expanded edition, University of Chicago Press, 1980.
• Allan Birnbaum, Some Latent Trait Models and Their Use in Inferring an Examinee's Ability, in Lord and Novick, Statistical Theories of Mental Test Scores, Addison-Wesley, 1968.
• R. D. Bock, M. Aitkin, Marginal Maximum Likelihood Estimation of Item Parameters: Application of an EM Algorithm, Psychometrika 46(4), 1981, pp. 443-459. paper
• J. Neyman, E. L. Scott, Consistent Estimates Based on Partially Consistent Observations, Econometrica 16(1), 1948. paper
• R. Philip Chalmers, mirt: A Multidimensional Item Response Theory Package for the R Environment, J. Stat. Softw. 48(6), 2012. paper · code
• J. Whitehill, T. Wu, J. Bergsma, J. Movellan, P. Ruvolo, Whose Vote Should Count More: Optimal Integration of Labels from Labelers of Unknown Expertise, NeurIPS 2009. paper
• B. Settles, G. T. LaFlair, M. Hagiwara, Machine Learning-Driven Language Assessment, TACL 8, 2020, pp. 247-263. paper
• J. P. Lalor, H. Wu, H. Yu, Building an Evaluation Scale using Item Response Theory, EMNLP 2016, pp. 648-657. paper
• HealthMeasures, Computer Adaptive Tests (CATs). page