From 73602fcfb778052c62ec5d5858258b29c7310846 Mon Sep 17 00:00:00 2001
From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com>
Date: Sat, 11 Jul 2026 09:23:29 +0100
Subject: [PATCH 1/2] test: dedicated chi-square correctness review + bug fixes
 (chi_square_test, chi_square_goodness_of_fit)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

PR #41 merged chi_square_test, chi_square_goodness_of_fit, and frequency_table
with only smoke-test coverage. This is the dedicated correctness review.

Bugs found and fixed in src/stats/inferential.jl:

1. Both functions computed the upper-tail p-value as `1 - cdf(...)`, which
   underflows to exactly 0.0 for large chi-squared statistics (verified
   interactively: cdf saturates to 1.0 in Float64 before the true tail
   probability reaches zero). Switched to `ccdf(...)`, accurate into the
   1e-22 range and beyond.
2. chi_square_test had no guard for r<2 or c<2: Chisq(0) throws an uncaught
   DomainError (verified interactively), and Cramer's V divides by zero via
   min(r,c)-1. Now a clean ArgumentError.
3. chi_square_test silently treated a fully zero row/column as if it didn't
   exist (via an `expected > 0` skip) while still charging its degree of
   freedom to df — an undefined comparison presented as a normal result.
   R's own chisq.test hits the same input and leaks NaN/NA; this now
   returns `nothing` + a "note" instead of fabricating a number.
4. chi_square_goodness_of_fit had no guard for k<2 categories (same
   Chisq(0) crash), no length/sum-to-1 validation on a caller-supplied
   expected_proportions, and no guard for a caller-supplied zero
   proportion ((O-0)^2/0 is Inf or NaN).
5. Neither function warned when expected cell counts fall below 5
   (Cochran's rule of thumb) — both now set a non-fatal "warning" key.
6. Added an optional, off-by-default Yates' continuity correction for 2x2
   tables (yates_correction=true), verified against R's
   chisq.test(correct=TRUE).

Also added require_nonnegative to src/stats/validation.jl (shared helper,
same pattern as require_positive).

frequency_table (src/stats/descriptive.jl) was reviewed and found
mathematically correct as-is — no production change needed there.

NEW test/chi_square_validation_test.jl: ground truth cross-checked against
both R 4.5.0's chisq.test and SciPy 1.15.3's chi2_contingency/chisquare
(agree to >=12 significant figures), covering an RxC independence example,
a non-uniform goodness-of-fit example, Yates correction, every degenerate
guard above (ArgumentError sites and nothing+note sites), the low-expected
warning, JSON-serialisable finiteness on normal cases, and executor-dispatch
parity. Wired into test/runtests.jl (include appended at the end of the
list per the task's overlap note).

Local verification: full suite green (flock-serialized WSL Julia run) —
424 (Full) + 155 (E2E) + 3800 (Property) + 62 (Reference) + 943 (Degenerate)
+ 80 (Guardrail) + 131 (Executor Router) + 208 (new Chi-Square Validation)
= 5803 passing, 0 failing. (One earlier run hit a pre-existing, unrelated
MANOVA test flake from an unseeded randn() call — confirmed unrelated by a
clean rerun with no code changes; not touched by this PR.)

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/stats/inferential.jl           | 134 +++++++++++--
 src/stats/validation.jl            |  10 +
 test/chi_square_validation_test.jl | 289 +++++++++++++++++++++++++++++
 test/runtests.jl                   |   4 +
 4 files changed, 418 insertions(+), 19 deletions(-)
 create mode 100644 test/chi_square_validation_test.jl

diff --git a/src/stats/inferential.jl b/src/stats/inferential.jl
index f2c9fb0..e2ec759 100644
--- a/src/stats/inferential.jl
+++ b/src/stats/inferential.jl
@@ -141,29 +141,77 @@ function one_way_anova(groups::Vector{Vector{Float64}}; alpha::Float64=0.05)
 end
 
 """
-    chi_square_test(observed::Matrix{Int}; alpha=0.05) -> Dict
+    chi_square_test(observed::Matrix{Int}; alpha=0.05, yates_correction=false) -> Dict
 
 CHI-SQUARE TEST OF INDEPENDENCE: Tests whether two categorical variables are
 independent in an r×c contingency table of observed frequencies.
-- EFFECT SIZE: Reports Cramér's V.
-- OUTPUT: chi-squared statistic, df = (r-1)(c-1), and p-value from the
-  chi-squared distribution.
+- EXPECTED COUNTS: `E_ij = row_sum_i * col_sum_j / n`.
+- STATISTIC: Pearson `χ² = Σ (O_ij - E_ij)² / E_ij`, `df = (r-1)(c-1)`.
+- P-VALUE: upper-tail `ccdf(Chisq(df), χ²)` (numerically stable for large χ²,
+  unlike `1 - cdf(...)` which underflows to exactly 0.0 in the far tail).
+- EFFECT SIZE: Reports Cramér's V = `√(χ² / (n·(min(r,c)-1)))`.
+- YATES' CONTINUITY CORRECTION: optional (`yates_correction=true`), valid only
+  for 2×2 tables: `χ² = Σ (|O-E| - 0.5)² / E`. Off by default to preserve the
+  standard (uncorrected) Pearson statistic; matches R's `chisq.test(correct=TRUE)`.
+- DEGENERATE GUARDS: throws `ArgumentError` for a malformed table (fewer than
+  2 rows/columns, negative counts, zero total observations — none of these
+  describe a valid r×c contingency table). Returns `nothing`+`"note"` (never
+  NaN/Inf) when the table shape is valid but a row or column sums to zero,
+  since the independence test is statistically undefined there (this is the
+  behaviour R's own `chisq.test` silently leaks as NaN/NA for the same input).
+  Emits a non-fatal `"warning"` when any expected count is below 5 (Cochran's
+  rule of thumb — the χ² approximation may be unreliable, but the statistic
+  is still returned).
 """
-function chi_square_test(observed::Matrix{Int}; alpha::Float64=0.05)
+function chi_square_test(observed::Matrix{Int}; alpha::Float64=0.05,
+                          yates_correction::Bool=false)
     r, c = size(observed)
+    require_at_least(r, 2, "observed rows")
+    require_at_least(c, 2, "observed columns")
+    require_nonnegative(observed, "observed")
+    if yates_correction && (r != 2 || c != 2)
+        throw(ArgumentError("yates_correction is only defined for 2×2 tables, got $(r)×$(c)"))
+    end
+
     n = sum(observed)
-    row_sums = sum(observed, dims=2)
-    col_sums = sum(observed, dims=1)
-
-    chi2 = 0.0
-    for i in 1:r, j in 1:c
-        expected = row_sums[i] * col_sums[j] / n
-        if expected > 0
-            chi2 += (observed[i, j] - expected)^2 / expected
-        end
+    require_at_least(n, 1, "sum(observed)")
+
+    row_sums = vec(sum(observed, dims=2))
+    col_sums = vec(sum(observed, dims=1))
+
+    # DEGENERATE GUARD: a row or column that sums to zero means that category
+    # was never observed at all. The expected-count formula still evaluates
+    # to a clean 0 for every cell in it (0 * col_sum / n = 0), but the
+    # resulting χ² would silently omit that category's contribution while
+    # still charging its degree of freedom — an undefined comparison, not a
+    # legitimate result. R's `chisq.test` hits the same condition and leaks
+    # `NaN`/`NA`; we return `nothing` + a note instead of fabricating a number.
+    if any(==(0), row_sums) || any(==(0), col_sums)
+        return Dict{String,Any}(
+            "chi_squared" => nothing, "df" => nothing, "p_value" => nothing,
+            "significant" => false, "cramers_v" => nothing,
+            "n" => n, "test_type" => "Chi-square test of independence",
+            "note" => "test undefined: a row or column has zero total observations",
+            "warning" => nothing
+        )
+    end
+
+    expected = [row_sums[i] * col_sums[j] / n for i in 1:r, j in 1:c]
+
+    chi2 = if yates_correction
+        sum((abs(observed[i, j] - expected[i, j]) - 0.5)^2 / expected[i, j]
+            for i in 1:r, j in 1:c)
+    else
+        sum((observed[i, j] - expected[i, j])^2 / expected[i, j] for i in 1:r, j in 1:c)
     end
     df = (r - 1) * (c - 1)
-    p_value = 1 - cdf(Chisq(df), chi2)
+    p_value = ccdf(Chisq(df), chi2)
+
+    min_expected = minimum(expected)
+    n_low = count(<(5), expected)
+    warning = min_expected < 5 ?
+        "$n_low of $(r*c) expected cell counts are below 5; the χ² approximation may be unreliable (Cochran's rule of thumb)" :
+        nothing
 
     return Dict{String,Any}(
         "chi_squared" => chi2,
@@ -172,7 +220,11 @@ function chi_square_test(observed::Matrix{Int}; alpha::Float64=0.05)
         "significant" => p_value < alpha,
         "cramers_v" => sqrt(chi2 / (n * (min(r, c) - 1))),
         "n" => n,
-        "test_type" => "Chi-square test of independence"
+        "test_type" => yates_correction ?
+            "Chi-square test of independence (Yates-corrected)" :
+            "Chi-square test of independence",
+        "note" => nothing,
+        "warning" => warning
     )
 end
 
@@ -182,18 +234,60 @@ end
 CHI-SQUARE GOODNESS-OF-FIT: Tests whether observed category counts match an
 expected distribution. Defaults to a uniform distribution across categories
 when `expected_proportions` is not supplied.
+- STATISTIC: `χ² = Σ (O_i - E_i)² / E_i` where `E_i = n * p_i`, `df = k - 1`.
+- P-VALUE: upper-tail `ccdf(Chisq(df), χ²)` (see `chi_square_test` docstring).
+- DEGENERATE GUARDS: throws `ArgumentError` for fewer than 2 categories,
+  negative counts, `expected_proportions` of mismatched length or not
+  summing to 1 (±1e-9), or zero total observations. Returns `nothing`+`"note"`
+  (never NaN/Inf) if any expected count is exactly zero (only reachable via a
+  caller-supplied zero proportion — the default uniform split cannot produce
+  this). Emits a non-fatal `"warning"` when any expected count is below 5.
 """
 function chi_square_goodness_of_fit(observed::Vector{Int},
                                     expected_proportions::Union{Vector{Float64},Nothing}=nothing;
                                     alpha::Float64=0.05)
     k = length(observed)
+    require_at_least(k, 2, "observed categories")
+    require_nonnegative(observed, "observed")
+
+    props = if isnothing(expected_proportions)
+        fill(1.0 / k, k)
+    else
+        require_equal_length(observed, expected_proportions, "observed", "expected_proportions")
+        require_nonnegative(expected_proportions, "expected_proportions")
+        isapprox(sum(expected_proportions), 1.0; atol=1e-9) || throw(ArgumentError(
+            "expected_proportions must sum to 1, got $(sum(expected_proportions))"))
+        expected_proportions
+    end
+
     n = sum(observed)
-    props = isnothing(expected_proportions) ? fill(1.0 / k, k) : expected_proportions
+    require_at_least(n, 1, "sum(observed)")
+
     expected = n .* props
 
+    # DEGENERATE GUARD: only reachable when a caller supplies a 0.0
+    # proportion for some category. (observed - 0)^2 / 0 is Inf (if that
+    # category was nonetheless observed) or NaN (0/0, if it wasn't) —
+    # either way, an undefined comparison, not a legitimate statistic.
+    if any(==(0.0), expected)
+        return Dict{String,Any}(
+            "chi_squared" => nothing, "df" => nothing, "p_value" => nothing,
+            "significant" => false, "expected" => expected,
+            "n" => n, "test_type" => "Chi-square goodness-of-fit",
+            "note" => "test undefined: expected_proportions assigns zero probability to an observed category",
+            "warning" => nothing
+        )
+    end
+
     chi2 = sum((observed .- expected) .^ 2 ./ expected)
     df = k - 1
-    p_value = 1 - cdf(Chisq(df), chi2)
+    p_value = ccdf(Chisq(df), chi2)
+
+    min_expected = minimum(expected)
+    n_low = count(<(5), expected)
+    warning = min_expected < 5 ?
+        "$n_low of $k expected cell counts are below 5; the χ² approximation may be unreliable (Cochran's rule of thumb)" :
+        nothing
 
     return Dict{String,Any}(
         "chi_squared" => chi2,
@@ -202,6 +296,8 @@ function chi_square_goodness_of_fit(observed::Vector{Int},
         "significant" => p_value < alpha,
         "expected" => expected,
         "n" => n,
-        "test_type" => "Chi-square goodness-of-fit"
+        "test_type" => "Chi-square goodness-of-fit",
+        "note" => nothing,
+        "warning" => warning
     )
 end
diff --git a/src/stats/validation.jl b/src/stats/validation.jl
index bd4a436..750784d 100644
--- a/src/stats/validation.jl
+++ b/src/stats/validation.jl
@@ -87,6 +87,16 @@ function require_probability(v::Real, name::String="value")
     return nothing
 end
 
+"""
+    require_nonnegative(x, name::String)
+
+Throw `ArgumentError` unless every element of `x` is `>= 0`.
+"""
+function require_nonnegative(x, name::String)
+    all(v -> v >= 0, x) || throw(ArgumentError("$name must not contain negative values"))
+    return nothing
+end
+
 """
     require_square(M::AbstractMatrix, name::String="matrix")
 
diff --git a/test/chi_square_validation_test.jl b/test/chi_square_validation_test.jl
new file mode 100644
index 0000000..75af7e2
--- /dev/null
+++ b/test/chi_square_validation_test.jl
@@ -0,0 +1,289 @@
+# SPDX-License-Identifier: MPL-2.0
+# Copyright (c) 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
+#
+# chi_square_validation_test.jl — Dedicated correctness review of the three
+# functions PR #41 merged with only smoke-test coverage: `chi_square_test`
+# and `chi_square_goodness_of_fit` (src/stats/inferential.jl), and the
+# completed `frequency_table` (src/stats/descriptive.jl). All three produce
+# numbers end users see verbatim, so this file checks the mathematics
+# against ground truth computed by a SECOND, independent implementation
+# rather than the library's own output.
+#
+# Ground truth was derived locally (2026-07-11, WSL Debian) with:
+#   - R 4.5.0's `chisq.test` (stats package, base R)
+#   - SciPy 1.15.3's `scipy.stats.chi2_contingency` / `chisquare`
+# Both implementations agreed to >= 12 significant figures on every case
+# below; only the constants (not the derivation scripts) are committed.
+#
+# BUGS FOUND AND FIXED during this review (see src/stats/inferential.jl):
+#   1. `chi_square_test` / `chi_square_goodness_of_fit` used `1 - cdf(...)`
+#      for the upper-tail p-value, which underflows to exactly 0.0 for
+#      large chi-squared statistics (verified: cdf saturates to 1.0 in
+#      Float64 well before the true tail probability reaches zero).
+#      Switched to `ccdf(...)`, which stays numerically accurate into the
+#      1e-22 range and beyond.
+#   2. `chi_square_test` had no guard for r<2 or c<2 columns/rows or a
+#      zero-observation table: `Chisq(0)` throws an uncaught `DomainError`
+#      (verified interactively), and `Cramér's V`'s `min(r,c)-1` divides by
+#      zero for a 1-row/1-column table. Both are now `ArgumentError`s.
+#   3. `chi_square_test` silently treated any table with a fully zero row
+#      or column as if that category didn't exist (via an `expected > 0`
+#      skip), while still charging its degree of freedom to `df` — an
+#      undefined comparison presented as a normal result. R's own
+#      `chisq.test` hits the same input and leaks `NaN`/`NA`; this function
+#      now returns `nothing` + a `"note"` explaining why, never a number.
+#   4. `chi_square_goodness_of_fit` had no guard for k<2 categories (same
+#      `Chisq(0)` crash as above), no length/sum-to-1 validation on a
+#      caller-supplied `expected_proportions`, and no guard for a
+#      caller-supplied zero proportion (`(O-0)^2/0` is `Inf` or `NaN`).
+#   5. Neither function warned when expected cell counts fall below 5
+#      (Cochran's rule of thumb) — both now set a non-fatal `"warning"`.
+#   6. Added an optional, explicitly-off-by-default Yates' continuity
+#      correction for 2x2 tables (`yates_correction=true`), verified
+#      against R's `chisq.test(..., correct=TRUE)`.
+#
+# `frequency_table` (src/stats/descriptive.jl) was reviewed and found
+# mathematically correct as-is: category counts, relative/cumulative
+# frequencies, and mode selection all check out against a hand count.
+# No production code change was needed there; it's still covered below.
+
+using Test
+using Statistikles
+
+@testset "Chi-Square Dedicated Validation" begin
+
+    # ═══════════════════════════════════════════════════════════════════
+    # 1. CHI-SQUARE TEST OF INDEPENDENCE — 3x3 ground truth
+    # ═══════════════════════════════════════════════════════════════════
+    @testset "Independence (3x3) vs R chisq.test / SciPy chi2_contingency" begin
+        # obs <- matrix(c(10,20,30, 6,9,17, 8,15,20), nrow=3, byrow=TRUE)
+        # R:     chisq.test(obs, correct=FALSE)
+        #        X-squared = 0.515205667618567, df = 4, p-value = 0.972003923728684
+        # SciPy: scipy.stats.chi2_contingency(obs, correction=False)
+        #        chi2 = 0.5152056676185671, dof = 4, p = 0.9720039237286836
+        # (R and SciPy agree to 12 significant figures; SciPy's extra ULP is
+        # float64 rounding noise, well inside the tolerance below.)
+        observed = [10 20 30; 6 9 17; 8 15 20]
+        r = Statistikles.chi_square_test(observed)
+
+        @test isapprox(r["chi_squared"], 0.515205667618567; atol = 1e-9)
+        @test r["df"] == 4
+        @test isapprox(r["p_value"], 0.972003923728684; atol = 1e-9)
+        @test r["significant"] == false
+        @test r["n"] == 135
+
+        # Cramér's V = sqrt(chi2 / (n * (min(r,c)-1))) = sqrt(0.515205667618567 / (135*2))
+        expected_cramers_v = sqrt(0.515205667618567 / (135 * 2))
+        @test isapprox(r["cramers_v"], expected_cramers_v; atol = 1e-9)
+
+        @test r["note"] === nothing
+        @test r["warning"] === nothing  # all expected counts are well above 5
+    end
+
+    @testset "Independence (2x2, uncorrected) vs R / SciPy" begin
+        # obs2 <- matrix(c(8,12,15,5), nrow=2, byrow=TRUE)
+        # R:     chisq.test(obs2, correct=FALSE)
+        #        X-squared = 5.012787723785166, df = 1, p-value = 0.0251607592004087
+        # SciPy: chi2_contingency(obs2, correction=False) agrees to 12 s.f.
+        observed = [8 12; 15 5]
+        r = Statistikles.chi_square_test(observed)
+        @test isapprox(r["chi_squared"], 5.012787723785166; atol = 1e-9)
+        @test r["df"] == 1
+        @test isapprox(r["p_value"], 0.0251607592004087; atol = 1e-9)
+        @test r["significant"] == true  # p < 0.05
+    end
+
+    @testset "Independence (2x2, Yates-corrected) vs R chisq.test(correct=TRUE)" begin
+        # R: chisq.test(obs2, correct=TRUE)
+        #    X-squared = 3.682864450127878, df = 1, p-value = 0.0549743287216943
+        observed = [8 12; 15 5]
+        r = Statistikles.chi_square_test(observed; yates_correction = true)
+        @test isapprox(r["chi_squared"], 3.682864450127878; atol = 1e-9)
+        @test r["df"] == 1
+        @test isapprox(r["p_value"], 0.0549743287216943; atol = 1e-9)
+        @test r["significant"] == false  # p > 0.05 (correction pulls it back from significance)
+        @test occursin("Yates", r["test_type"])
+    end
+
+    @testset "yates_correction rejected on a non-2x2 table" begin
+        observed = [10 20 30; 6 9 17; 8 15 20]
+        @test_throws ArgumentError Statistikles.chi_square_test(observed; yates_correction = true)
+    end
+
+    # ═══════════════════════════════════════════════════════════════════
+    # 2. CHI-SQUARE GOODNESS-OF-FIT — non-uniform expected proportions
+    # ═══════════════════════════════════════════════════════════════════
+    @testset "Goodness-of-fit (non-uniform) vs R chisq.test(p=...) / SciPy chisquare" begin
+        # observed_gof <- c(18, 22, 15, 25, 20); props <- c(.15,.25,.10,.30,.20)
+        # R:     chisq.test(observed_gof, p=props)
+        #        X-squared = 4.293333333333333, df = 4, p-value = 0.367760644428913
+        # SciPy: scipy.stats.chisquare(observed_gof, f_exp=n*props)
+        #        statistic = 4.293333333333333, pvalue = 0.36776064442891276
+        observed = [18, 22, 15, 25, 20]
+        props = [0.15, 0.25, 0.10, 0.30, 0.20]
+        r = Statistikles.chi_square_goodness_of_fit(observed, props)
+
+        @test isapprox(r["chi_squared"], 4.293333333333333; atol = 1e-9)
+        @test r["df"] == 4
+        @test isapprox(r["p_value"], 0.367760644428913; atol = 1e-9)
+        @test r["significant"] == false
+        @test isapprox(r["expected"], [15.0, 25.0, 10.0, 30.0, 20.0]; atol = 1e-12)
+        @test r["n"] == 100
+        @test r["note"] === nothing
+        @test r["warning"] === nothing  # min expected = 10, above the Cochran threshold
+    end
+
+    @testset "Goodness-of-fit (uniform default) closed-form cross-check" begin
+        # k=3 uniform categories, n=30 -> expected=10 each. Deviations
+        # [+5,0,-5] -> chi2 = (25+0+25)/10 = 5.0 exactly, df=2.
+        # df=2 has a closed-form chi-square survival function independent
+        # of the library's own Distributions.cdf call path:
+        #   P(chi^2_2 > x) = exp(-x/2)
+        observed = [15, 10, 5]
+        r = Statistikles.chi_square_goodness_of_fit(observed)
+        @test r["chi_squared"] == 5.0
+        @test r["df"] == 2
+        @test isapprox(r["p_value"], exp(-2.5); atol = 1e-12)
+    end
+
+    # ═══════════════════════════════════════════════════════════════════
+    # 3. DEGENERATE-INPUT HANDLING — clean errors, no NaN/Inf leaks
+    # ═══════════════════════════════════════════════════════════════════
+    @testset "chi_square_test: 1-row table throws ArgumentError (was: DomainError)" begin
+        @test_throws ArgumentError Statistikles.chi_square_test(reshape([10, 20, 30], 1, 3))
+    end
+
+    @testset "chi_square_test: negative counts throw ArgumentError" begin
+        @test_throws ArgumentError Statistikles.chi_square_test([10 -5; 3 8])
+    end
+
+    @testset "chi_square_test: all-zero table throws ArgumentError" begin
+        @test_throws ArgumentError Statistikles.chi_square_test([0 0; 0 0])
+    end
+
+    @testset "chi_square_test: zero row total -> nothing + note, not NaN" begin
+        observed = [0 0 0; 6 9 17; 8 15 20]
+        r = Statistikles.chi_square_test(observed)
+        @test r["chi_squared"] === nothing
+        @test r["df"] === nothing
+        @test r["p_value"] === nothing
+        @test r["cramers_v"] === nothing
+        @test r["significant"] == false
+        @test r["note"] !== nothing
+        @test occursin("zero total", r["note"])
+        assert_finite_and_serialisable(r)
+    end
+
+    @testset "chi_square_test: zero column total -> nothing + note, not NaN" begin
+        observed = [0 20 30; 0 9 17; 0 15 20]
+        r = Statistikles.chi_square_test(observed)
+        @test r["chi_squared"] === nothing
+        @test r["note"] !== nothing
+        assert_finite_and_serialisable(r)
+    end
+
+    @testset "chi_square_goodness_of_fit: k=1 category throws ArgumentError" begin
+        @test_throws ArgumentError Statistikles.chi_square_goodness_of_fit([42])
+    end
+
+    @testset "chi_square_goodness_of_fit: mismatched proportions length throws ArgumentError" begin
+        @test_throws ArgumentError Statistikles.chi_square_goodness_of_fit([10, 20, 30], [0.5, 0.5])
+    end
+
+    @testset "chi_square_goodness_of_fit: proportions not summing to 1 throws ArgumentError" begin
+        @test_throws ArgumentError Statistikles.chi_square_goodness_of_fit([10, 20, 30], [0.5, 0.5, 0.5])
+    end
+
+    @testset "chi_square_goodness_of_fit: negative counts throw ArgumentError" begin
+        @test_throws ArgumentError Statistikles.chi_square_goodness_of_fit([10, -5, 30])
+    end
+
+    @testset "chi_square_goodness_of_fit: zero-observation table throws ArgumentError" begin
+        @test_throws ArgumentError Statistikles.chi_square_goodness_of_fit([0, 0, 0])
+    end
+
+    @testset "chi_square_goodness_of_fit: zero expected proportion -> nothing + note, not Inf" begin
+        # Category 3 is assigned zero probability but was observed anyway:
+        # (O-E)^2/E = (5-0)^2/0 = Inf under the naive formula.
+        observed = [10, 15, 5]
+        props = [0.5, 0.5, 0.0]
+        r = Statistikles.chi_square_goodness_of_fit(observed, props)
+        @test r["chi_squared"] === nothing
+        @test r["note"] !== nothing
+        assert_finite_and_serialisable(r)
+    end
+
+    @testset "chi_square_test: low expected counts trigger a warning, not an error" begin
+        # Small counts spread over a 3x2 table so several expected cells
+        # fall below 5, but the table itself is still well-formed.
+        observed = [2 1; 1 2; 3 1]
+        r = Statistikles.chi_square_test(observed)
+        @test r["chi_squared"] !== nothing  # statistic is still returned
+        @test r["warning"] !== nothing
+        @test occursin("below 5", r["warning"])
+        assert_finite_and_serialisable(r)
+    end
+
+    # ═══════════════════════════════════════════════════════════════════
+    # 4. JSON-SERIALISABLE FINITENESS on a normal case (both functions)
+    # ═══════════════════════════════════════════════════════════════════
+    @testset "chi_square_test: normal case is fully finite + JSON-serialisable" begin
+        r = Statistikles.chi_square_test([10 20 30; 6 9 17; 8 15 20])
+        assert_finite_and_serialisable(r)
+    end
+
+    @testset "chi_square_goodness_of_fit: normal case is fully finite + JSON-serialisable" begin
+        r = Statistikles.chi_square_goodness_of_fit([18, 22, 15, 25, 20],
+                                                      [0.15, 0.25, 0.10, 0.30, 0.20])
+        assert_finite_and_serialisable(r)
+    end
+
+    # ═══════════════════════════════════════════════════════════════════
+    # 5. FREQUENCY_TABLE — reviewed, found correct; covered here too
+    # ═══════════════════════════════════════════════════════════════════
+    @testset "frequency_table: hand-counted ground truth" begin
+        data = ["x", "y", "x", "x", "z", "y", "x"]
+        # Hand count: x×4, y×2, z×1, n=7; categories sort alphabetically.
+        r = Statistikles.frequency_table(data)
+        @test r["categories"] == ["x", "y", "z"]
+        @test r["frequencies"] == [4, 2, 1]
+        @test isapprox(r["relative_frequencies"], [400 / 7, 200 / 7, 100 / 7]; atol = 1e-9)
+        @test r["cumulative_frequencies"] == [4, 6, 7]
+        @test isapprox(r["cumulative_relative_frequencies"], [400 / 7, 600 / 7, 100.0]; atol = 1e-9)
+        @test r["n"] == 7
+        @test r["n_categories"] == 3
+        @test r["mode"] == "x"
+        assert_finite_and_serialisable(r)
+    end
+
+    @testset "frequency_table: empty input returns a clean error, not a crash" begin
+        r = Statistikles.frequency_table(String[])
+        @test haskey(r, "error")
+        assert_finite_and_serialisable(r)
+    end
+
+    # ═══════════════════════════════════════════════════════════════════
+    # 6. EXECUTOR DISPATCH — chi_square tool still matches the direct call
+    #    (guards against the fix above breaking the router's Matrix{Int}
+    #    coercion path in src/tools/executor.jl)
+    # ═══════════════════════════════════════════════════════════════════
+    @testset "Executor dispatch: chi_square independence unaffected by the guards" begin
+        observed_rows = [[10, 20, 30], [6, 9, 17], [8, 15, 20]]
+        direct = Statistikles.chi_square_test([10 20 30; 6 9 17; 8 15 20])
+        via_tool = Statistikles.execute_tool("chi_square",
+            Dict{String,Any}("type" => "independence", "observed" => observed_rows))
+        @test isapprox(via_tool["chi_squared"], direct["chi_squared"]; atol = 1e-12)
+        @test isapprox(via_tool["p_value"], direct["p_value"]; atol = 1e-12)
+    end
+
+    @testset "Executor dispatch: chi_square goodness-of-fit unaffected by the guards" begin
+        direct = Statistikles.chi_square_goodness_of_fit([18, 22, 15, 25, 20],
+                                                           [0.15, 0.25, 0.10, 0.30, 0.20])
+        via_tool = Statistikles.execute_tool("chi_square",
+            Dict{String,Any}("type" => "goodness_of_fit", "observed" => [18, 22, 15, 25, 20],
+                              "expected_proportions" => [0.15, 0.25, 0.10, 0.30, 0.20]))
+        @test isapprox(via_tool["chi_squared"], direct["chi_squared"]; atol = 1e-12)
+        @test isapprox(via_tool["p_value"], direct["p_value"]; atol = 1e-12)
+    end
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index c3ede8e..c169f0a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1353,3 +1353,7 @@ include("guardrail_test.jl")
 # Executor router coverage: every LLM-facing tool name in definitions.jl
 # is exercised through execute_tool (or explicitly skipped with a reason)
 include("executor_router_test.jl")
+
+# Chi-square dedicated correctness review: chi_square_test, chi_square_goodness_of_fit,
+# and frequency_table checked against independently-computed ground truth (R / SciPy)
+include("chi_square_validation_test.jl")

From 61c31a5539f30284dc2f2fda48b9fb0d440293bc Mon Sep 17 00:00:00 2001
From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com>
Date: Sat, 11 Jul 2026 09:36:26 +0100
Subject: [PATCH 2/2] fix: clamp Yates continuity correction per-cell to fix
 chi2 inflation

Reviewer found (must_fix): chi_square_test's yates_correction branch
(src/stats/inferential.jl) computed each cell's correction as
(abs(O-E) - 0.5)^2 / E with no max(0.0, ...) clamp. Whenever every cell's
|O-E| < 0.5 (near-independence 2x2 tables), this squares a *negative*
number into a *positive* correction, inflating chi2 instead of pulling
it toward zero -- the opposite of what Yates' correction is for.

Verified failure: observed=[[10,10],[10,11]] (common |O-E|=0.2439<0.5).
- R 4.5.0 chisq.test(correct=TRUE):            X-squared = 0, p = 1
- SciPy 1.15.3 chi2_contingency(correction=True): chi2 = 0.0,  p = 1.0
- Buggy formula (pre-fix):                      chi2 = 0.0256, p != 1

Fixed to max(0.0, abs(O-E) - 0.5)^2 / E (equivalently: subtract
min(0.5, |O-E|) before squaring), matching R and SciPy exactly.

Added a regression test in test/chi_square_validation_test.jl for this
exact case, since the existing Yates test ([[8,12],[15,5]], |O-E|=3.5)
has every cell's |O-E| > 0.5 and cannot distinguish the correct formula
from the buggy one -- both agree there.

Local verification: full suite green (flock-serialized WSL Julia run) --
424 (Full) + 155 (E2E) + 3800 (Property) + 62 (Reference) + 943
(Degenerate) + 80 (Guardrail) + 131 (Executor Router) + 213 (Chi-Square
Validation, was 208) = 5808 passing, 0 failing.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 src/stats/inferential.jl           | 15 ++++++++---
 test/chi_square_validation_test.jl | 40 ++++++++++++++++++++++++++++++
 2 files changed, 52 insertions(+), 3 deletions(-)

diff --git a/src/stats/inferential.jl b/src/stats/inferential.jl
index e2ec759..9f369ad 100644
--- a/src/stats/inferential.jl
+++ b/src/stats/inferential.jl
@@ -151,8 +151,12 @@ independent in an r×c contingency table of observed frequencies.
   unlike `1 - cdf(...)` which underflows to exactly 0.0 in the far tail).
 - EFFECT SIZE: Reports Cramér's V = `√(χ² / (n·(min(r,c)-1)))`.
 - YATES' CONTINUITY CORRECTION: optional (`yates_correction=true`), valid only
-  for 2×2 tables: `χ² = Σ (|O-E| - 0.5)² / E`. Off by default to preserve the
-  standard (uncorrected) Pearson statistic; matches R's `chisq.test(correct=TRUE)`.
+  for 2×2 tables: `χ² = Σ (max(0, |O-E| - 0.5))² / E`. The `max(0, ...)` clamp
+  matters whenever every cell's `|O-E| < 0.5` (near-independence tables): without
+  it the unclamped formula can *inflate* χ² instead of correcting it downward.
+  Off by default to preserve the standard (uncorrected) Pearson statistic;
+  matches R's `chisq.test(correct=TRUE)` and SciPy's
+  `chi2_contingency(correction=True)`.
 - DEGENERATE GUARDS: throws `ArgumentError` for a malformed table (fewer than
   2 rows/columns, negative counts, zero total observations — none of these
   describe a valid r×c contingency table). Returns `nothing`+`"note"` (never
@@ -199,7 +203,12 @@ function chi_square_test(observed::Matrix{Int}; alpha::Float64=0.05,
     expected = [row_sums[i] * col_sums[j] / n for i in 1:r, j in 1:c]
 
     chi2 = if yates_correction
-        sum((abs(observed[i, j] - expected[i, j]) - 0.5)^2 / expected[i, j]
+        # Clamp each cell's correction to a minimum of 0: without
+        # `max(0.0, ...)`, a cell with |O-E| < 0.5 would square a *negative*
+        # number, inflating χ² instead of correcting it downward (e.g.
+        # observed=[10 10; 10 11] has |O-E|=0.2439 in every cell and must
+        # yield χ²=0.0, not 0.0256 — see test/chi_square_validation_test.jl).
+        sum((max(0.0, abs(observed[i, j] - expected[i, j]) - 0.5))^2 / expected[i, j]
             for i in 1:r, j in 1:c)
     else
         sum((observed[i, j] - expected[i, j])^2 / expected[i, j] for i in 1:r, j in 1:c)
diff --git a/test/chi_square_validation_test.jl b/test/chi_square_validation_test.jl
index 75af7e2..262718e 100644
--- a/test/chi_square_validation_test.jl
+++ b/test/chi_square_validation_test.jl
@@ -41,6 +41,17 @@
 #   6. Added an optional, explicitly-off-by-default Yates' continuity
 #      correction for 2x2 tables (`yates_correction=true`), verified
 #      against R's `chisq.test(..., correct=TRUE)`.
+#   7. [Adversarial-review follow-up] The initial Yates implementation
+#      computed `(abs(O-E) - 0.5)^2 / E` with no `max(0.0, ...)` clamp, so
+#      any cell with `|O-E| < 0.5` squared a *negative* correction into a
+#      *positive* one, inflating chi2 instead of correcting it downward —
+#      confirmed to disagree with both R's `chisq.test(correct=TRUE)` and
+#      SciPy's `chi2_contingency(correction=True)` on
+#      `observed=[[10,10],[10,11]]` (buggy chi2 ≈ 0.0256 vs ground truth
+#      0.0). Fixed to `max(0.0, abs(O-E) - 0.5)^2 / E`; regression test
+#      added below ("near-independence" Yates case) since the original
+#      Yates test only used `[[8,12],[15,5]]` where `|O-E|=3.5`, which
+#      cannot distinguish the correct formula from the buggy one.
 #
 # `frequency_table` (src/stats/descriptive.jl) was reviewed and found
 # mathematically correct as-is: category counts, relative/cumulative
@@ -105,6 +116,35 @@ using Statistikles
         @test occursin("Yates", r["test_type"])
     end
 
+    @testset "Independence (2x2, Yates-corrected, near-independence) vs R / SciPy — regression for missing clamp" begin
+        # This case specifically guards the Yates continuity-correction clamp:
+        # every cell's |O-E| = 0.2439 < 0.5, so the *correct* formula
+        # max(0, |O-E| - 0.5)^2 / E must clamp every term to zero, giving
+        # chi2 = 0.0 exactly. An earlier version of chi_square_test omitted
+        # the clamp (computed (|O-E| - 0.5)^2 / E directly), which squares a
+        # small negative number into a *positive* contribution instead of
+        # zero, inflating chi2 to ~0.0256 and deflating the p-value — enough
+        # to flip a significance decision near alpha in general, though not
+        # for this particular alpha=0.05 case. The previous Yates test above
+        # ([[8,12],[15,5]], |O-E|=3.5) does NOT exercise the clamp, since
+        # both the correct and buggy formulas agree once |O-E| > 0.5.
+        #
+        # obs3 <- matrix(c(10,10,10,11), nrow=2, byrow=TRUE)
+        # R 4.5.0:     chisq.test(obs3, correct=TRUE)
+        #              X-squared = 0, df = 1, p-value = 1
+        # SciPy 1.15.3: chi2_contingency(obs3, correction=True)
+        #              chi2 = 0.0, p = 1.0, dof = 1
+        #              expected = [[9.75609756, 10.24390244],
+        #                          [10.24390244, 10.75609756]]
+        observed = [10 10; 10 11]
+        r = Statistikles.chi_square_test(observed; yates_correction = true)
+        @test isapprox(r["chi_squared"], 0.0; atol = 1e-9)
+        @test r["df"] == 1
+        @test isapprox(r["p_value"], 1.0; atol = 1e-9)
+        @test r["significant"] == false
+        @test occursin("Yates", r["test_type"])
+    end
+
     @testset "yates_correction rejected on a non-2x2 table" begin
         observed = [10 20 30; 6 9 17; 8 15 20]
         @test_throws ArgumentError Statistikles.chi_square_test(observed; yates_correction = true)
