Skip to content

Commit a040d48

Browse files
committed
jc: scale-invariant tolerances, exact omega regressions, independent sign oracle
Addresses all three findings from the #889 review (4856837957) and the #888 post-merge note (5181899915). Each was reproduced before fixing. P1 — omega was still UNIT-DEPENDENT. Four absolute floors survived in omega_total: `.max(1.0)` in the lambda^2 guard, the psi guard, and both sign cutoffs. Below unit variance those degrade to an absolute 1e-9, so on data whose covariances are ~1e-16 every tolerance is enormous, every sign comparison vacuous, and (measured) omega flips 0.25 -> 0.75 purely by rescaling. All four removed; tolerances are now relative to the item's own variance, and covariance cutoffs to sqrt(var_i * var_j). This is the SAME defect class as the R-squared absolute pivot fixed two commits earlier. I fixed the instance and left the class — the floors were in the function I was actively editing. Constant items (zero variance) are now rejected explicitly: they carry no signal, leave the model unidentified, and would make every scale-relative tolerance degenerate. P1 — the regressions blessed the old wrong answer. The zero-loading test asserted only `(0.0..=1.0).contains(&w)`, which 0.75 also satisfies, so it proved "did not reject" rather than "computed the signed omega". Worse, the invariance tests compared against a SELF-COMPUTED baseline, so a consistently wrong implementation would satisfy them at every scale and permutation. Both absolute values are now pinned as constants derived by hand: ZERO_ANCHOR_OMEGA = 0.25 (lambda = [0,+1,-1,+1], psi = [V,0,V,V]) CONGENERIC_OMEGA = 18/19 (exact rational; the 8-digit decimal is 1.05e-9 off and a 1e-9 invariance tolerance correctly rejects it) and the rescale/permutation tests assert against those, not against themselves. Invariance alone never pins a value. P2 — the "independent" frustration oracle mirrored the implementation. It reconstructed signs from each anchor row and looked for conflicts, which is an anchor-shaped decision procedure, not the specification. Replaced with the literal predicate: fix the free global sign and enumerate all 2^(k-1) = 8 assignments, asserting NONE satisfies every materially non-zero covariance edge, plus an anti-vacuity assert that the sign graph carries >= 4 material edges. Also removed (review): the global-observation-sign-flip test. Negating every item leaves the covariance matrix bit-identical, and omega_total consumes only that matrix, so the assertion could not fail for ANY implementation — the workspace's own "assertion implied by the code it tests". Replaced by common rescaling, item permutation (including rotations that move the zero-loading item into index 0), and a single-item flip that must CHANGE omega. 118 lib + 13 doctests green; stats.rs clippy-clean. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01K3RyLEbuNSHxxB3NTTrGki
1 parent 5602bdc commit a040d48

1 file changed

Lines changed: 170 additions & 53 deletions

File tree

‎crates/jc/src/stats.rs‎

Lines changed: 170 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -386,6 +386,13 @@ pub fn omega_total(items: &[Vec<f64>]) -> Option<f64> {
386386
}
387387
}
388388

389+
// A constant item has zero variance: it carries no signal, makes every
390+
// scale-relative tolerance below degenerate, and leaves the model
391+
// unidentified. Reject rather than absorb it.
392+
if (0..k).any(|i| cov[i][i] <= 0.0) {
393+
return None;
394+
}
395+
389396
// λ_i² averaged over every admissible triad (j,k ≠ i, j ≠ k, σ_jk ≠ 0).
390397
let mut lambda = vec![0.0f64; k];
391398
for i in 0..k {
@@ -418,7 +425,13 @@ pub fn omega_total(items: &[Vec<f64>]) -> Option<f64> {
418425
// has λ² == σ_ii exactly in real arithmetic but lands a few ulps either
419426
// side of it in f64, and a bare `< 0.0` test would reject a valid model
420427
// as misfit. Only a violation LARGER than rounding is a real one.
421-
let tol = 1e-9 * cov[i][i].abs().max(1.0);
428+
// RELATIVE to the item's own variance, with NO absolute floor. A
429+
// `.max(1.0)` here is a statement about UNITS: on data whose
430+
// covariances are ~1e-16 it makes a 1e-9 tolerance enormous, every
431+
// comparison vacuous, and (measured) flipped ω from 0.25 back to
432+
// 0.75 purely by rescaling — the SAME defect class as the R²
433+
// absolute pivot, fixed there first and left here until review.
434+
let tol = 1e-9 * cov[i][i].abs();
422435
if lam_sq < -tol || !lam_sq.is_finite() {
423436
return None; // single-factor model violated (negative common variance)
424437
}
@@ -462,9 +475,7 @@ pub fn omega_total(items: &[Vec<f64>]) -> Option<f64> {
462475
// A materially-zero covariance means the anchor's or item j's loading
463476
// is ~0; a ~0-loading item contributes ~nothing to Σλ, so either sign
464477
// is harmless. Keep +1.
465-
let scale = (cov[anchor][anchor].abs() * cov[j][j].abs())
466-
.sqrt()
467-
.max(1.0);
478+
let scale = (cov[anchor][anchor].abs() * cov[j][j].abs()).sqrt();
468479
if cov[anchor][j] < -1e-9 * scale {
469480
sign[j] = -1.0;
470481
}
@@ -477,7 +488,7 @@ pub fn omega_total(items: &[Vec<f64>]) -> Option<f64> {
477488
// constraints; see the function's doc for what is and is not checked.)
478489
for i in 0..k {
479490
for j in (i + 1)..k {
480-
let scale = (cov[i][i].abs() * cov[j][j].abs()).sqrt().max(1.0);
491+
let scale = (cov[i][i].abs() * cov[j][j].abs()).sqrt();
481492
let cutoff = 1e-9 * scale;
482493
if cov[i][j].abs() <= cutoff {
483494
continue; // too near zero to carry a sign
@@ -496,7 +507,7 @@ pub fn omega_total(items: &[Vec<f64>]) -> Option<f64> {
496507
let mut sum_psi = 0.0;
497508
for i in 0..k {
498509
let psi = cov[i][i] - lambda[i] * lambda[i];
499-
let tol = 1e-9 * cov[i][i].abs().max(1.0);
510+
let tol = 1e-9 * cov[i][i].abs();
500511
if psi < -tol {
501512
// Heywood case: estimated common variance genuinely exceeds the
502513
// item's total variance → model misfit, not a reliability value.
@@ -1316,46 +1327,124 @@ mod tests {
13161327
);
13171328
}
13181329

1319-
#[test]
1320-
fn omega_is_invariant_to_a_global_sign_flip() {
1321-
// The one-factor model is identified only up to a GLOBAL flip, and ω
1322-
// uses (Σλ)², so flipping every item must leave ω unchanged. This also
1323-
// proves the sign anchor (s_0 = +1) introduces no arbitrary bias.
1324-
let items = congeneric_items();
1325-
let flipped: Vec<Vec<f64>> = items
1326-
.iter()
1327-
.map(|it| it.iter().map(|v| -v).collect())
1328-
.collect();
1329-
let a = omega_total(&items).unwrap();
1330-
let b = omega_total(&flipped).unwrap();
1331-
assert!(approx(a, b, 1e-12), "global flip changed ω: {a} vs {b}");
1330+
/// ω of [`congeneric_items`], EXACT: `(Σλ)² = 48`, `Σψ = 8/3`, so
1331+
/// `ω = 48 / (48 + 8/3) = 144/152 = 18/19`. Written as the rational rather
1332+
/// than a truncated decimal — an 8-digit literal is 1.05e-9 off, which a
1333+
/// 1e-9 invariance tolerance correctly rejects.
1334+
const CONGENERIC_OMEGA: f64 = 18.0 / 19.0;
1335+
1336+
/// ω of [`zero_anchor_items`], derived exactly below.
1337+
///
1338+
/// With f, e₀, eₐ, e_b mutually orthogonal ±1 vectors of length 8 (each of
1339+
/// sample variance V = 8/7):
1340+
/// loadings λ = [0, +1, −1, +1]·√V → Σλ = √V
1341+
/// common variance (Σλ)² = V
1342+
/// residuals ψ = [V, 0, V, V] → Σψ = 3V
1343+
/// ω = V / (V + 3V) = 0.25
1344+
/// The pre-fix UNSIGNED implementation returns 0.75 on this fixture, so
1345+
/// pinning 0.25 is what distinguishes a correct sign fix from "did not
1346+
/// reject" — a `(0.0..=1.0)` range check blesses both.
1347+
const ZERO_ANCHOR_OMEGA: f64 = 0.25;
1348+
1349+
/// λ = [0, +1, −1, +1]: item 0 carries NO loading, so it is the worst
1350+
/// possible sign anchor — the fixture that exposes anchor bias.
1351+
fn zero_anchor_items() -> Vec<Vec<f64>> {
1352+
let f = [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0];
1353+
let e0 = [1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0];
1354+
let ea = [1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0];
1355+
let eb = [1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0];
1356+
vec![
1357+
e0.to_vec(), // λ = 0
1358+
f.to_vec(), // λ = +1
1359+
(0..8).map(|i| -f[i] + ea[i]).collect::<Vec<_>>(), // λ = −1
1360+
(0..8).map(|i| f[i] + eb[i]).collect::<Vec<_>>(), // λ = +1
1361+
]
13321362
}
13331363

13341364
#[test]
1335-
fn omega_sign_anchor_survives_a_zero_loading_item() {
1336-
// REGRESSION: signs are read from ONE item's covariance row, so the
1337-
// choice of that row matters. Anchoring on item 0 read every sign off
1338-
// an all-zero row when item 0 carried no loading — every sign defaulted
1339-
// to +1, and the consistency check then manufactured conflicts and
1340-
// REJECTED this valid one-factor dataset (measured: 2 false conflicts).
1365+
fn omega_is_invariant_to_a_common_rescale() {
1366+
// REPLACES a vacuous "global sign flip" test (external review).
1367+
// Negating EVERY item leaves the covariance matrix bit-identical —
1368+
// `cov(−x,−y) = cov(x,y)` — and `omega_total` consumes only that
1369+
// matrix, so the old assertion could not fail for any implementation.
1370+
// It was the workspace's own "assertion implied by the code it tests".
13411371
//
1342-
// λ = [0, +1, +1, −1] with orthogonal residuals: item 0 is pure noise,
1343-
// the rest genuinely carry mixed signs. Σλ = 0+1+1−1 = 1, so the
1344-
// common variance is small but the model is real and must not be
1345-
// rejected.
1346-
let f = [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0];
1347-
let e0 = [1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0];
1348-
let e2 = [1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0];
1349-
let e3 = [1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0];
1350-
let items = vec![
1351-
e0.to_vec(), // λ = 0
1352-
f.to_vec(), // λ = +1
1353-
(0..8).map(|i| f[i] + e2[i]).collect::<Vec<f64>>(), // λ = +1
1354-
(0..8).map(|i| -f[i] + e3[i]).collect::<Vec<f64>>(), // λ = −1
1355-
];
1356-
let w = omega_total(&items)
1357-
.expect("a valid one-factor set must not be rejected over a weak anchor");
1358-
assert!((0.0..=1.0).contains(&w), "ω out of range: {w}");
1372+
// The real invariance: scaling every item by one constant sends
1373+
// λ → cλ and ψ → c²ψ, leaving the ratio unchanged. This DOES vary the
1374+
// input, and it is what the absolute `.max(1.0)` tolerance floors
1375+
// broke — measured, ω flipped 0.25 → 0.75 at scale 1e-8.
1376+
//
1377+
// Each fixture is checked against its KNOWN value, not against a
1378+
// self-computed baseline (external review P1): an implementation that
1379+
// is consistently WRONG across every scale would satisfy a
1380+
// self-referential invariance check. Invariance alone never pins a
1381+
// value.
1382+
for (fixture, expected) in [
1383+
(congeneric_items(), CONGENERIC_OMEGA),
1384+
(zero_anchor_items(), ZERO_ANCHOR_OMEGA),
1385+
] {
1386+
let reference = omega_total(&fixture).expect("baseline must fit");
1387+
assert!(
1388+
approx(reference, expected, 1e-9),
1389+
"baseline ω={reference}, expected {expected}"
1390+
);
1391+
for c in [1e-8, 1e-4, 1e-2, 1e2, 1e6] {
1392+
let scaled: Vec<Vec<f64>> = fixture
1393+
.iter()
1394+
.map(|it| it.iter().map(|v| v * c).collect())
1395+
.collect();
1396+
let got = omega_total(&scaled)
1397+
.unwrap_or_else(|| panic!("scale {c:e} rejected a valid model"));
1398+
assert!(
1399+
approx(got, reference, 1e-9),
1400+
"scale {c:e}: ω={got} vs {reference}"
1401+
);
1402+
}
1403+
}
1404+
}
1405+
1406+
#[test]
1407+
fn omega_is_invariant_to_item_permutation() {
1408+
// ω is a property of the item SET, not of their order — but signs are
1409+
// read from one anchor row, so a bad anchor rule shows up here. The
1410+
// rotation that puts the ZERO-loading item at index 0 is exactly the
1411+
// case that broke the previous anchor-on-item-0 implementation.
1412+
let base = zero_anchor_items();
1413+
// Pin the absolute value here too — see the rescale test's note.
1414+
let reference = ZERO_ANCHOR_OMEGA;
1415+
assert!(approx(omega_total(&base).unwrap(), reference, 1e-9));
1416+
let k = base.len();
1417+
for shift in 1..k {
1418+
let rotated: Vec<Vec<f64>> = (0..k).map(|i| base[(i + shift) % k].clone()).collect();
1419+
let got = omega_total(&rotated)
1420+
.unwrap_or_else(|| panic!("rotation by {shift} rejected a valid model"));
1421+
assert!(
1422+
approx(got, reference, 1e-9),
1423+
"shift {shift}: ω={got} vs {reference}"
1424+
);
1425+
}
1426+
// Reverse order too — a different permutation class from rotation.
1427+
let mut reversed = base.clone();
1428+
reversed.reverse();
1429+
assert!(approx(omega_total(&reversed).unwrap(), reference, 1e-9));
1430+
}
1431+
1432+
#[test]
1433+
fn omega_changes_when_one_item_is_flipped() {
1434+
// The non-vacuous half of the sign story: flipping ONE item genuinely
1435+
// changes the model (that item's loading changes sign), so ω MUST
1436+
// move. A sign-blind implementation would return the same value —
1437+
// which is precisely the bug this fixture family was built for.
1438+
let base = congeneric_items();
1439+
let reference = omega_total(&base).unwrap();
1440+
let mut flipped = base.clone();
1441+
flipped[1] = flipped[1].iter().map(|v| -v).collect();
1442+
let got =
1443+
omega_total(&flipped).expect("flipping one item is still a valid one-factor model");
1444+
assert!(
1445+
(got - reference).abs() > 0.05,
1446+
"one-item flip must move ω: {got} vs {reference}"
1447+
);
13591448
}
13601449

13611450
#[test]
@@ -1418,21 +1507,49 @@ mod tests {
14181507
"item {i}: ψ would trip the HEYWOOD guard"
14191508
);
14201509
}
1421-
// …and that the pattern is frustrated from EVERY anchor, so no anchor
1422-
// choice could rescue it.
1423-
for anchor in 0..k {
1510+
// …and that the pattern is genuinely FRUSTRATED, by the LITERAL
1511+
// specification rather than by the implementation's own procedure
1512+
// (external review P2). "Reconstruct signs from each anchor row and
1513+
// look for a conflict" is an implementation-shaped shortcut and is not
1514+
// the same predicate on sparse sign graphs. The specification is:
1515+
// does ANY assignment of item signs satisfy every materially non-zero
1516+
// covariance edge? At k = 4 that is 2^(k−1) = 8 candidates once the
1517+
// free global flip is fixed — small enough to enumerate exhaustively.
1518+
let material = |i: usize, j: usize| -> bool {
1519+
c[i][j].abs() > 1e-9 * (c[i][i].abs() * c[j][j].abs()).sqrt()
1520+
};
1521+
let edges: Vec<(usize, usize)> = (0..k)
1522+
.flat_map(|i| ((i + 1)..k).map(move |j| (i, j)))
1523+
.filter(|&(i, j)| material(i, j))
1524+
.collect();
1525+
// Anti-vacuity: frustration over an empty or near-empty edge set is
1526+
// meaningless, so the sign graph must actually carry constraints.
1527+
assert!(
1528+
edges.len() >= 4,
1529+
"only {} material edges — frustration would be vacuous",
1530+
edges.len()
1531+
);
1532+
let mut satisfying: Option<Vec<f64>> = None;
1533+
for mask in 0..(1u32 << (k - 1)) {
14241534
let mut sign = vec![1.0f64; k];
1425-
for j in 0..k {
1426-
if j != anchor && c[anchor][j] < 0.0 {
1427-
sign[j] = -1.0;
1535+
for (j, sj) in sign.iter_mut().enumerate().skip(1) {
1536+
if (mask >> (j - 1)) & 1 == 1 {
1537+
*sj = -1.0;
14281538
}
14291539
}
1430-
let conflicts = (0..k)
1431-
.flat_map(|i| ((i + 1)..k).map(move |j| (i, j)))
1432-
.filter(|&(i, j)| c[i][j] != 0.0 && c[i][j].signum() != sign[i] * sign[j])
1433-
.count();
1434-
assert!(conflicts > 0, "anchor {anchor} accepts it — not frustrated");
1540+
if edges
1541+
.iter()
1542+
.all(|&(i, j)| c[i][j].signum() == sign[i] * sign[j])
1543+
{
1544+
satisfying = Some(sign);
1545+
break;
1546+
}
14351547
}
1548+
assert!(
1549+
satisfying.is_none(),
1550+
"fixture is NOT frustrated — {satisfying:?} satisfies all {} edges",
1551+
edges.len()
1552+
);
14361553

14371554
// Therefore only the sign-consistency guard can produce the rejection.
14381555
assert_eq!(omega_total(&items), None);

0 commit comments

Comments
 (0)