Add Thurstone (1927) Case V paired-comparison scaling

CHANGELOG.md

@@ -129,6 +129,19 @@
 
 ### Added
 
+- **Thurstone Case V paired-comparison scaling (Thurstone, 1927, as
+  implemented by psych's `thurstone()`; Revelle, 2025 — source READ)**
+  (`fast_mlsirm.thurstone_case_v`; in Rust
+  `mlsirm_core::scaling::thurstone_case_v`): scales n objects from an n×n
+  choice-probability matrix (`choice[i, j]` = P(column j preferred over row
+  i)) via `scale_j = colmean(qnorm(choice))_j − min`, fitted model
+  `Phi(scale_j − scale_i)`, residuals, and psych's goodness of fit
+  `1 − sse/ssc` over the full model matrix (pinning the psych *code*
+  behavior; the .Rd "lower off diagonal" prose is stale). Entries must be
+  strictly in (0, 1) — a deliberate safety divergence from psych, whose
+  direct path admits infinite normal quantiles. Pinned against a 50-digit
+  mpmath oracle on asymmetric, exactly-consistent, and intransitive
+  fixtures.
 - **Composite linking (Holland & Strawderman, 2011; as cited by Albano,
   2016, JSS 74(8), eqs. 31-32)** (`fast_mlsirm.composite_linking`; in Rust
   `mlsirm_core::equating::composite_linking`): weighted average of H

crates/fast-mlsirm-py/src/lib.rs

@@ -5121,6 +5121,26 @@ fn composite_linking(
     Ok(d.into())
 }
 
+/// Thurstone (1927) Case V paired-comparison scaling, as implemented by
+/// psych's `thurstone()` (see `mlsirm_core::scaling`). `choice` is a flat
+/// row-major n*n matrix; returns dict with scale, gof, model, residual.
+#[pyfunction]
+#[pyo3(signature = (choice, n))]
+fn thurstone_case_v(
+    py: Python<'_>,
+    choice: PyReadonlyArray1<'_, f64>,
+    n: usize,
+) -> PyResult<Py<pyo3::types::PyDict>> {
+    let res = mlsirm_core::scaling::thurstone_case_v(choice.as_slice()?, n)
+        .map_err(PyValueError::new_err)?;
+    let d = pyo3::types::PyDict::new(py);
+    d.set_item("scale", PyArray1::from_slice(py, &res.scale))?;
+    d.set_item("gof", res.gof)?;
+    d.set_item("model", PyArray1::from_slice(py, &res.model))?;
+    d.set_item("residual", PyArray1::from_slice(py, &res.residual))?;
+    Ok(d.into())
+}
+
 /// GPCM/nominal softmax cell log-probabilities at one node (parity surface for
 /// the NumPy `category_logprobs` reference).
 #[pyfunction]
@@ -7457,6 +7477,7 @@ fn fast_mlsirm_core(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(circle_arc_equate, m)?)?;
     m.add_function(wrap_pyfunction!(nominal_weights_mean_equate, m)?)?;
     m.add_function(wrap_pyfunction!(composite_linking, m)?)?;
+    m.add_function(wrap_pyfunction!(thurstone_case_v, m)?)?;
     m.add_function(wrap_pyfunction!(circle_arc_middle_anchor, m)?)?;
     m.add_function(wrap_pyfunction!(loglinear_smooth, m)?)?;
     m.add_function(wrap_pyfunction!(person_fit_stat, m)?)?;

crates/mlsirm-core/src/lib.rs

@@ -34,6 +34,7 @@ pub mod reliability;
 pub mod rsm;
 pub mod rt;
 pub mod rt_joint;
+pub mod scaling;
 pub mod scoring;
 pub mod security;
 pub mod standard_setting;

crates/mlsirm-core/src/scaling.rs

@@ -0,0 +1,157 @@
+//! Thurstone Case V paired-comparison scaling.
+//!
+//! Implements Thurstone's (1927) law of comparative judgment under the
+//! Case V assumptions (equal and uncorrelated discriminal dispersions)
+//! exactly as implemented by the CRAN `psych` package's `thurstone()`
+//! main routine.
+//!
+//! Source status:
+//! - **READ**: `psych` (Revelle) `R/thurstone.R` (main routine) and
+//!   `R/matrix.addition.R` (the `%+%` operator used to form the model
+//!   matrix), fetched from CRAN. Every formula below is traceable to
+//!   those lines.
+//! - **NOT READ (as-cited)**: Thurstone, L. L. (1927). A law of
+//!   comparative judgment. *Psychological Review, 34*, 273-286.
+//!   https://doi.org/10.1037/h0070288 — paywalled at verification time.
+//!   Attribution is therefore "Thurstone (1927), Case V, as implemented
+//!   in psych's `thurstone()`" (same governance precedent as the
+//!   Guttman (1945) module in `reliability.rs`).
+//!
+//! Algorithm (thurstone.R lines 25-28), with `choice[i][j]` = proportion
+//! of judgments preferring column object `j` over row object `i`
+//! (psych convention, `man/thurstone.Rd`):
+//!
+//! ```text
+//! z_ij     = qnorm(choice_ij)
+//! S_j      = mean_i(z_ij) - min_k mean_i(z_ik)   # column means, min scale = 0
+//! model_ij = pnorm(S_j - S_i)
+//! error    = model - choice
+//! GF       = 1 - sum(error^2) / sum(choice^2)    # FULL matrix incl. diagonal
+//! ```
+//!
+//! DERIVED (verified by reading `matrix.addition.R` lines 11-17): in
+//! `pnorm(-scale.values %+% t(scale.values))` the `%+%` operator coerces
+//! `-scale.values` to an n-by-1 column whose row sums are `-S_i` and
+//! `t(scale.values)` is 1-by-n with column sums `S_j`, so the argument is
+//! `S_j - S_i`.
+//!
+//! DOC/CODE CONFLICT (documented): `man/thurstone.Rd` describes a
+//! "lower off diagonal" goodness of fit, but the code returns a single GF
+//! over the FULL matrix with no `lower.tri` masking. This implementation
+//! pins CODE behavior.
+//!
+//! DELIBERATE SAFETY DIVERGENCE: psych's direct-matrix path passes 0/1
+//! proportions straight into `qnorm`, producing infinities; this
+//! implementation rejects entries outside the open interval (0, 1)
+//! instead (psych itself clamps to [1/(4n), 1-1/(4n)] in its upstream
+//! `choice.mat()` converter, which is out of scope here along with the
+//! rank-order/item converters).
+//!
+//! Numerics: reuses the shared approximate kernels
+//! [`crate::nodes::inv_normal_cdf`] (Acklam, relative error < 1.15e-9)
+//! and [`crate::fitstats::erfc`] (Numerical Recipes rational,
+//! |error| < 1.2e-7). Outputs inherit that approximation error; the test
+//! pins are verified to agree with a 50-digit oracle within 1e-6 on the
+//! tested fixtures (no universal accuracy claim, especially for
+//! proportions near 0 or 1).
+//!
+//! # References
+//!
+//! Thurstone, L. L. (1927). A law of comparative judgment. *Psychological
+//! Review, 34*, 273-286. https://doi.org/10.1037/h0070288
+//!
+//! Revelle, W. (2024). *psych: Procedures for psychological, psychometric,
+//! and personality research* (R package). Northwestern University.
+//! https://CRAN.R-project.org/package=psych
+
+use crate::fitstats::erfc;
+use crate::nodes::inv_normal_cdf;
+
+/// Result of a Thurstone Case V scaling.
+#[derive(Clone, Debug)]
+pub struct ThurstoneResult {
+    /// Scale values, length `n`, shifted so the minimum is exactly 0
+    /// (thurstone.R line 25).
+    pub scale: Vec<f64>,
+    /// Goodness of fit `1 - sum(residual^2)/sum(choice^2)` over the full
+    /// matrix including the diagonal (thurstone.R line 28).
+    pub gof: f64,
+    /// Model choice matrix `pnorm(S_j - S_i)`, row-major `n * n`.
+    pub model: Vec<f64>,
+    /// Residual matrix `model - choice`, row-major `n * n`.
+    pub residual: Vec<f64>,
+}
+
+fn norm_cdf(z: f64) -> f64 {
+    0.5 * erfc(-z / std::f64::consts::SQRT_2)
+}
+
+/// Thurstone Case V scaling of a square choice matrix.
+///
+/// `choice` is row-major `n * n`; `choice[i*n + j]` is the proportion of
+/// judgments preferring column object `j` over row object `i`. Entries
+/// must lie strictly in (0, 1); the diagonal is used as-is (psych does
+/// not special-case it).
+pub fn thurstone_case_v(choice: &[f64], n: usize) -> Result<ThurstoneResult, String> {
+    if n < 2 {
+        return Err("thurstone_case_v needs at least 2 objects".into());
+    }
+    if choice.len() != n * n {
+        return Err(format!(
+            "choice must be a row-major {n}x{n} matrix ({} entries), got {}",
+            n * n,
+            choice.len()
+        ));
+    }
+    for &c in choice {
+        if !c.is_finite() {
+            return Err("choice proportions must be finite".into());
+        }
+        if c <= 0.0 || c >= 1.0 {
+            return Err("choice proportions must lie strictly in (0, 1)".into());
+        }
+    }
+
+    // Column means of z = qnorm(choice), shifted so min scale value is 0
+    // (thurstone.R line 25).
+    let mut colmeans = vec![0.0f64; n];
+    for i in 0..n {
+        for j in 0..n {
+            colmeans[j] += inv_normal_cdf(choice[i * n + j]);
+        }
+    }
+    for m in colmeans.iter_mut() {
+        *m /= n as f64;
+    }
+    let min = colmeans.iter().cloned().fold(f64::INFINITY, f64::min);
+    let scale: Vec<f64> = colmeans.iter().map(|&m| m - min).collect();
+
+    // model_ij = pnorm(S_j - S_i) (thurstone.R line 26 via %+%, DERIVED
+    // above); residual and full-matrix GF (lines 27-28).
+    let mut model = vec![0.0f64; n * n];
+    let mut residual = vec![0.0f64; n * n];
+    let mut sse = 0.0f64;
+    let mut ssc = 0.0f64;
+    for i in 0..n {
+        for j in 0..n {
+            let m = norm_cdf(scale[j] - scale[i]);
+            let e = m - choice[i * n + j];
+            model[i * n + j] = m;
+            residual[i * n + j] = e;
+            sse += e * e;
+            ssc += choice[i * n + j] * choice[i * n + j];
+        }
+    }
+    let gof = 1.0 - sse / ssc;
+
+    Ok(ThurstoneResult {
+        scale,
+        gof,
+        model,
+        residual,
+    })
+}
+
+#[cfg(test)]
+#[path = "../../../tests/unit/scaling_tests.rs"]
+mod tests;

python/fast_mlsirm/__init__.py

@@ -22,6 +22,7 @@
 from .linking import irt_link as irt_link, IrtLinkResult as IrtLinkResult
 from .equating import equate_observed_scores as equate_observed_scores, equate_neat as equate_neat, EquateResult as EquateResult, equate_observed_scores_kernel as equate_observed_scores_kernel, loglinear_smooth as loglinear_smooth, equate_neat_linear as equate_neat_linear, equating_standard_errors as equating_standard_errors, circle_arc_equate as circle_arc_equate, circle_arc_middle_anchor as circle_arc_middle_anchor, CircleArcResult as CircleArcResult, nominal_weights_mean_equate as nominal_weights_mean_equate, composite_linking as composite_linking
 from .rt import fit_response_times as fit_response_times, RtFit as RtFit, fit_speed_accuracy as fit_speed_accuracy, rt_person_fit as rt_person_fit
+from .scaling import thurstone_case_v as thurstone_case_v, ThurstoneResult as ThurstoneResult
 from .cdm import fit_cdm as fit_cdm, CdmFit as CdmFit, fit_gdina as fit_gdina, GdinaFit as GdinaFit, validate_q_matrix as validate_q_matrix, QMatrixValidation as QMatrixValidation, gdina_wald_selection as gdina_wald_selection, WaldModelSelection as WaldModelSelection, fit_ho_cdm as fit_ho_cdm, HoCdmFit as HoCdmFit, fit_ho_gdina as fit_ho_gdina, HoGdinaFit as HoGdinaFit, fit_seq_gdina as fit_seq_gdina, SeqGdinaFit as SeqGdinaFit, fit_seq_gdina_qr as fit_seq_gdina_qr, SeqGdinaQrFit as SeqGdinaQrFit
 from .mixture import fit_mixture as fit_mixture, MixtureFit as MixtureFit
 from .crm import fit_crm as fit_crm, CrmFit as CrmFit
@@ -225,6 +226,8 @@
     "circle_arc_equate",
     "nominal_weights_mean_equate",
     "composite_linking",
+    "thurstone_case_v",
+    "ThurstoneResult",
     "circle_arc_middle_anchor",
     "CircleArcResult",
     "fit_response_times",

python/fast_mlsirm/scaling.py

@@ -0,0 +1,72 @@
+"""Thurstone (1927) Case V paired-comparison scaling.
+
+Thin wrapper over the Rust core (``mlsirm_core::scaling``). The algorithm
+follows the ``thurstone()`` function of the psych R package (Revelle, 2025;
+source READ), which implements Thurstone's (1927) law of comparative judgment
+under Case V (equal discriminal dispersions, zero correlations). Thurstone
+(1927) itself was NOT read; it is cited as the origin of the model as
+described by the psych source.
+
+References
+----------
+Revelle, W. (2025). psych: Procedures for psychological, psychometric, and
+    personality research (R package). https://CRAN.R-project.org/package=psych
+Thurstone, L. L. (1927). A law of comparative judgment. Psychological
+    Review, 34(4), 273-286. [NOT READ; cited as described in psych source]
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+
+@dataclass
+class ThurstoneResult:
+    """Case V scaling result: ``scale[j]`` is the scale value of object j
+    (minimum shifted to exactly 0), ``gof`` the psych goodness of fit
+    ``1 - sse/ssc`` over the FULL model matrix (including the diagonal --
+    this pins the psych *code* behavior; the .Rd prose saying "lower off
+    diagonal" is stale), ``model`` the fitted choice probabilities
+    ``Phi(scale[j] - scale[i])`` and ``residual = model - choice``, both
+    shaped (n, n)."""
+
+    scale: np.ndarray
+    gof: float
+    model: np.ndarray
+    residual: np.ndarray
+
+
+def thurstone_case_v(choice) -> ThurstoneResult:
+    """Scale n objects from an n x n choice-probability matrix.
+
+    ``choice[i, j]`` is the proportion of judges preferring object *j* over
+    object *i* (psych convention: column beats row). All entries must be
+    strictly in (0, 1) -- a deliberate safety divergence from psych, whose
+    direct path lets ``qnorm(0)/qnorm(1)`` produce infinities.
+    """
+    from .fitstats import _core_module
+
+    arr = np.asarray(choice)
+    if np.iscomplexobj(arr):
+        raise ValueError("thurstone_case_v: choice must be real-valued")
+    if arr.dtype == object:
+        try:
+            arr = arr.astype(np.float64)
+        except (TypeError, ValueError) as exc:
+            raise ValueError(
+                "thurstone_case_v: choice must be numeric"
+            ) from exc
+    arr = np.ascontiguousarray(arr, dtype=np.float64)
+    if arr.ndim != 2 or arr.shape[0] != arr.shape[1]:
+        raise ValueError("thurstone_case_v: choice must be a square 2-D matrix")
+    n = arr.shape[0]
+    core = _core_module()
+    res = core.thurstone_case_v(arr.ravel(), n)
+    return ThurstoneResult(
+        scale=np.asarray(res["scale"], dtype=np.float64),
+        gof=float(res["gof"]),
+        model=np.asarray(res["model"], dtype=np.float64).reshape(n, n),
+        residual=np.asarray(res["residual"], dtype=np.float64).reshape(n, n),
+    )

tests/test_paper_features.py

@@ -8538,4 +8538,70 @@ def test_errors(self):
         with pytest.raises(ValueError):
             composite_linking([t + 1j], [1.0])
         with pytest.raises(ValueError):
-            composite_linking([np.array([1.0, "a"], dtype=object)], [1.0])
+            composite_linking([np.array([1.0, "a"], dtype=object)], [1.0])
+
+class TestThurstoneCaseV:
+    """Thurstone (1927) Case V scaling via psych's thurstone() algorithm.
+
+    Pins from files/thurstone_oracle.py (mpmath 50-digit, EXECUTED).
+    All asserts read the ThurstoneResult returned by the crate.
+    """
+
+    def test_fixture_a_scale_gof_model(self):
+        # Kills sign/colmean/min-shift mutants: scale + GF + off-diagonal
+        # model cells pinned against the 50-digit oracle (asymmetric,
+        # nonzero-residual fixture).
+        import numpy as np
+        from fast_mlsirm import thurstone_case_v
+
+        choice = np.array(
+            [[0.5, 0.7, 0.9], [0.3, 0.5, 0.8], [0.1, 0.2, 0.5]]
+        )
+        r = thurstone_case_v(choice)
+        np.testing.assert_allclose(
+            r.scale,
+            [0.0, 0.49624378579592261, 1.3097082924567186],
+            rtol=0,
+            atol=1e-6,
+        )
+        assert abs(r.gof - 0.99986967677023893) < 1e-6
+        assert abs(r.model[0, 2] - 0.9048527313906016) < 1e-6
+        assert abs(r.model[2, 0] - 0.095147268609398397) < 1e-6
+        assert abs(r.residual[0, 1] - (-0.00986121100602988)) < 1e-6
+        assert r.model.shape == (3, 3) and r.residual.shape == (3, 3)
+
+    def test_fixture_b_consistent_roundtrip(self):
+        # Exactly Case-V-consistent 2x2: GF must be 1, model reproduces
+        # the input (round-trip anchor; variance mutations are covered by
+        # fixture A/C pins, not this zero-residual fixture).
+        import numpy as np
+        from fast_mlsirm import thurstone_case_v
+
+        choice = np.array([[0.5, 0.75], [0.25, 0.5]])
+        r = thurstone_case_v(choice)
+        assert abs(r.gof - 1.0) < 1e-9
+        np.testing.assert_allclose(
+            r.scale, [0.0, 0.67448975019608174], rtol=0, atol=1e-6
+        )
+        np.testing.assert_allclose(r.model, choice, rtol=0, atol=1e-6)
+
+    def test_errors(self):
+        import numpy as np
+        import pytest
+        from fast_mlsirm import thurstone_case_v
+
+        good = np.array([[0.5, 0.75], [0.25, 0.5]])
+        with pytest.raises(ValueError):
+            thurstone_case_v(good[:1])  # non-square
+        with pytest.raises(ValueError):
+            thurstone_case_v(np.array([[0.5]]))  # n < 2
+        bad = good.copy()
+        bad[0, 1] = 1.0  # boundary excluded: strictly in (0, 1)
+        with pytest.raises(ValueError):
+            thurstone_case_v(bad)
+        with pytest.raises(ValueError):
+            thurstone_case_v(good + 0j)  # complex rejected
+        with pytest.raises(ValueError):
+            thurstone_case_v(
+                np.array([[0.5, "a"], [0.25, 0.5]], dtype=object)
+            )