Skip to content

Commit f488ce4

Browse files
test: dedicated chi-square correctness review + Yates-clamp fix (chi_square_test, chi_square_goodness_of_fit) (#47)
## What & why PR #41 merged `chi_square_test` and `chi_square_goodness_of_fit` (`src/stats/inferential.jl`), plus the completed `frequency_table` (`src/stats/descriptive.jl`), with only smoke-test coverage. These produce numbers end users see verbatim, so this is the dedicated correctness review the wave-1 audit called for (`.claude/tasks/prod-readiness/w2-5-reference-validation.md` lineage / split of PR #41's follow-up). ### Bugs found and fixed in `src/stats/inferential.jl` 1. **`1 - cdf(...)` underflow** — both functions computed the upper-tail p-value as `1 - cdf(...)`, which underflows to exactly `0.0` for large χ² 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`, and Cramér'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`/`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 2×2 tables (`yates_correction=true`). 7. **[This push — reviewer must_fix]** The Yates implementation from item 6 computed `(abs(O-E) - 0.5)^2 / E` with **no `max(0.0, ...)` clamp**. Whenever every cell's `|O-E| < 0.5` (near-independence 2×2 tables), this squares a *negative* number into a *positive* correction, **inflating** χ² instead of correcting it downward — the opposite of what Yates' correction is for. Verified: `observed=[[10,10],[10,11]]` (common `|O-E|=0.2439<0.5`) gave buggy χ²≈0.0256 vs R/SciPy ground truth `0.0`. Fixed to `max(0.0, abs(O-E) - 0.5)^2 / E` (equivalently subtract `min(0.5, |O-E|)` before squaring), matching R's `chisq.test(correct=TRUE)` and SciPy's `chi2_contingency(correction=True)` exactly. `frequency_table` (`src/stats/descriptive.jl`) was reviewed and found mathematically correct as-is — no production change needed there. ### Ground truth sources Derived locally (2026-07-11, WSL Debian) with: - **R 4.5.0**'s `chisq.test` (base `stats` package) - **SciPy 1.15.3**'s `scipy.stats.chi2_contingency` / `chisquare` Both agree to ≥12 significant figures on every case in the new test file; only the constants (not the derivation scripts) are committed, each with a source comment. New: the Yates-clamp regression case (`[[10,10],[10,11]]`) was independently re-verified against both R 4.5.0 and SciPy 1.15.3 in this session (`chisq.test(obs, correct=TRUE)` → `X-squared = 0, df = 1, p-value = 1`; `chi2_contingency(obs, correction=True)` → `chi2=0.0, p=1.0, dof=1`). ## Test coverage (`test/chi_square_validation_test.jl`, wired into `test/runtests.jl`) - RxC independence example (3×3) vs R/SciPy. - 2×2 independence, uncorrected, vs R/SciPy. - 2×2 independence, Yates-corrected (`[[8,12],[15,5]]`, `|O-E|=3.5`), vs R. - **New**: 2×2 independence, Yates-corrected, near-independence (`[[10,10],[10,11]]`, `|O-E|=0.2439<0.5` in every cell) — this is the case that catches the missing clamp; the existing `[[8,12],[15,5]]` case does not, since both the correct and buggy formulas agree once `|O-E| > 0.5`. - `yates_correction` rejected on non-2×2 tables. - Goodness-of-fit, non-uniform proportions, vs R/SciPy. - Every degenerate guard (`ArgumentError` sites and `nothing`+`"note"` sites), the low-expected-count warning, JSON-serialisable finiteness on normal cases, and executor-dispatch parity. ## Verification run Full suite, flock-serialized WSL Julia run (`julia --project=. -e 'using Pkg; Pkg.test()'`): ``` Statistikles Full Test Suite 424/424 Statistikles E2E Pipeline Tests 155/155 Statistikles Property-Based Tests 3800/3800 Reference Validation (ground truth) 62/62 Degenerate Input Guards 943/943 Neural-Boundary Guardrail 80/80 Executor Router Coverage 131/131 Chi-Square Dedicated Validation 213/213 (was 208 before this push's regression test) ----------------------------------------- Total: 5808 passing, 0 failing ``` ## Skipped / out of scope None. `frequency_table` review found no bug (documented above, not re-litigated here). 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
1 parent 49b8300 commit f488ce4

4 files changed

Lines changed: 467 additions & 19 deletions

File tree

src/stats/inferential.jl

Lines changed: 124 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -141,29 +141,86 @@ function one_way_anova(groups::Vector{Vector{Float64}}; alpha::Float64=0.05)
141141
end
142142

143143
"""
144-
chi_square_test(observed::Matrix{Int}; alpha=0.05) -> Dict
144+
chi_square_test(observed::Matrix{Int}; alpha=0.05, yates_correction=false) -> Dict
145145
146146
CHI-SQUARE TEST OF INDEPENDENCE: Tests whether two categorical variables are
147147
independent in an r×c contingency table of observed frequencies.
148-
- EFFECT SIZE: Reports Cramér's V.
149-
- OUTPUT: chi-squared statistic, df = (r-1)(c-1), and p-value from the
150-
chi-squared distribution.
148+
- EXPECTED COUNTS: `E_ij = row_sum_i * col_sum_j / n`.
149+
- STATISTIC: Pearson `χ² = Σ (O_ij - E_ij)² / E_ij`, `df = (r-1)(c-1)`.
150+
- P-VALUE: upper-tail `ccdf(Chisq(df), χ²)` (numerically stable for large χ²,
151+
unlike `1 - cdf(...)` which underflows to exactly 0.0 in the far tail).
152+
- EFFECT SIZE: Reports Cramér's V = `√(χ² / (n·(min(r,c)-1)))`.
153+
- YATES' CONTINUITY CORRECTION: optional (`yates_correction=true`), valid only
154+
for 2×2 tables: `χ² = Σ (max(0, |O-E| - 0.5))² / E`. The `max(0, ...)` clamp
155+
matters whenever every cell's `|O-E| < 0.5` (near-independence tables): without
156+
it the unclamped formula can *inflate* χ² instead of correcting it downward.
157+
Off by default to preserve the standard (uncorrected) Pearson statistic;
158+
matches R's `chisq.test(correct=TRUE)` and SciPy's
159+
`chi2_contingency(correction=True)`.
160+
- DEGENERATE GUARDS: throws `ArgumentError` for a malformed table (fewer than
161+
2 rows/columns, negative counts, zero total observations — none of these
162+
describe a valid r×c contingency table). Returns `nothing`+`"note"` (never
163+
NaN/Inf) when the table shape is valid but a row or column sums to zero,
164+
since the independence test is statistically undefined there (this is the
165+
behaviour R's own `chisq.test` silently leaks as NaN/NA for the same input).
166+
Emits a non-fatal `"warning"` when any expected count is below 5 (Cochran's
167+
rule of thumb — the χ² approximation may be unreliable, but the statistic
168+
is still returned).
151169
"""
152-
function chi_square_test(observed::Matrix{Int}; alpha::Float64=0.05)
170+
function chi_square_test(observed::Matrix{Int}; alpha::Float64=0.05,
171+
yates_correction::Bool=false)
153172
r, c = size(observed)
173+
require_at_least(r, 2, "observed rows")
174+
require_at_least(c, 2, "observed columns")
175+
require_nonnegative(observed, "observed")
176+
if yates_correction && (r != 2 || c != 2)
177+
throw(ArgumentError("yates_correction is only defined for 2×2 tables, got $(r)×$(c)"))
178+
end
179+
154180
n = sum(observed)
155-
row_sums = sum(observed, dims=2)
156-
col_sums = sum(observed, dims=1)
157-
158-
chi2 = 0.0
159-
for i in 1:r, j in 1:c
160-
expected = row_sums[i] * col_sums[j] / n
161-
if expected > 0
162-
chi2 += (observed[i, j] - expected)^2 / expected
163-
end
181+
require_at_least(n, 1, "sum(observed)")
182+
183+
row_sums = vec(sum(observed, dims=2))
184+
col_sums = vec(sum(observed, dims=1))
185+
186+
# DEGENERATE GUARD: a row or column that sums to zero means that category
187+
# was never observed at all. The expected-count formula still evaluates
188+
# to a clean 0 for every cell in it (0 * col_sum / n = 0), but the
189+
# resulting χ² would silently omit that category's contribution while
190+
# still charging its degree of freedom — an undefined comparison, not a
191+
# legitimate result. R's `chisq.test` hits the same condition and leaks
192+
# `NaN`/`NA`; we return `nothing` + a note instead of fabricating a number.
193+
if any(==(0), row_sums) || any(==(0), col_sums)
194+
return Dict{String,Any}(
195+
"chi_squared" => nothing, "df" => nothing, "p_value" => nothing,
196+
"significant" => false, "cramers_v" => nothing,
197+
"n" => n, "test_type" => "Chi-square test of independence",
198+
"note" => "test undefined: a row or column has zero total observations",
199+
"warning" => nothing
200+
)
201+
end
202+
203+
expected = [row_sums[i] * col_sums[j] / n for i in 1:r, j in 1:c]
204+
205+
chi2 = if yates_correction
206+
# Clamp each cell's correction to a minimum of 0: without
207+
# `max(0.0, ...)`, a cell with |O-E| < 0.5 would square a *negative*
208+
# number, inflating χ² instead of correcting it downward (e.g.
209+
# observed=[10 10; 10 11] has |O-E|=0.2439 in every cell and must
210+
# yield χ²=0.0, not 0.0256 — see test/chi_square_validation_test.jl).
211+
sum((max(0.0, abs(observed[i, j] - expected[i, j]) - 0.5))^2 / expected[i, j]
212+
for i in 1:r, j in 1:c)
213+
else
214+
sum((observed[i, j] - expected[i, j])^2 / expected[i, j] for i in 1:r, j in 1:c)
164215
end
165216
df = (r - 1) * (c - 1)
166-
p_value = 1 - cdf(Chisq(df), chi2)
217+
p_value = ccdf(Chisq(df), chi2)
218+
219+
min_expected = minimum(expected)
220+
n_low = count(<(5), expected)
221+
warning = min_expected < 5 ?
222+
"$n_low of $(r*c) expected cell counts are below 5; the χ² approximation may be unreliable (Cochran's rule of thumb)" :
223+
nothing
167224

168225
return Dict{String,Any}(
169226
"chi_squared" => chi2,
@@ -172,7 +229,11 @@ function chi_square_test(observed::Matrix{Int}; alpha::Float64=0.05)
172229
"significant" => p_value < alpha,
173230
"cramers_v" => sqrt(chi2 / (n * (min(r, c) - 1))),
174231
"n" => n,
175-
"test_type" => "Chi-square test of independence"
232+
"test_type" => yates_correction ?
233+
"Chi-square test of independence (Yates-corrected)" :
234+
"Chi-square test of independence",
235+
"note" => nothing,
236+
"warning" => warning
176237
)
177238
end
178239

@@ -182,18 +243,60 @@ end
182243
CHI-SQUARE GOODNESS-OF-FIT: Tests whether observed category counts match an
183244
expected distribution. Defaults to a uniform distribution across categories
184245
when `expected_proportions` is not supplied.
246+
- STATISTIC: `χ² = Σ (O_i - E_i)² / E_i` where `E_i = n * p_i`, `df = k - 1`.
247+
- P-VALUE: upper-tail `ccdf(Chisq(df), χ²)` (see `chi_square_test` docstring).
248+
- DEGENERATE GUARDS: throws `ArgumentError` for fewer than 2 categories,
249+
negative counts, `expected_proportions` of mismatched length or not
250+
summing to 1 (±1e-9), or zero total observations. Returns `nothing`+`"note"`
251+
(never NaN/Inf) if any expected count is exactly zero (only reachable via a
252+
caller-supplied zero proportion — the default uniform split cannot produce
253+
this). Emits a non-fatal `"warning"` when any expected count is below 5.
185254
"""
186255
function chi_square_goodness_of_fit(observed::Vector{Int},
187256
expected_proportions::Union{Vector{Float64},Nothing}=nothing;
188257
alpha::Float64=0.05)
189258
k = length(observed)
259+
require_at_least(k, 2, "observed categories")
260+
require_nonnegative(observed, "observed")
261+
262+
props = if isnothing(expected_proportions)
263+
fill(1.0 / k, k)
264+
else
265+
require_equal_length(observed, expected_proportions, "observed", "expected_proportions")
266+
require_nonnegative(expected_proportions, "expected_proportions")
267+
isapprox(sum(expected_proportions), 1.0; atol=1e-9) || throw(ArgumentError(
268+
"expected_proportions must sum to 1, got $(sum(expected_proportions))"))
269+
expected_proportions
270+
end
271+
190272
n = sum(observed)
191-
props = isnothing(expected_proportions) ? fill(1.0 / k, k) : expected_proportions
273+
require_at_least(n, 1, "sum(observed)")
274+
192275
expected = n .* props
193276

277+
# DEGENERATE GUARD: only reachable when a caller supplies a 0.0
278+
# proportion for some category. (observed - 0)^2 / 0 is Inf (if that
279+
# category was nonetheless observed) or NaN (0/0, if it wasn't) —
280+
# either way, an undefined comparison, not a legitimate statistic.
281+
if any(==(0.0), expected)
282+
return Dict{String,Any}(
283+
"chi_squared" => nothing, "df" => nothing, "p_value" => nothing,
284+
"significant" => false, "expected" => expected,
285+
"n" => n, "test_type" => "Chi-square goodness-of-fit",
286+
"note" => "test undefined: expected_proportions assigns zero probability to an observed category",
287+
"warning" => nothing
288+
)
289+
end
290+
194291
chi2 = sum((observed .- expected) .^ 2 ./ expected)
195292
df = k - 1
196-
p_value = 1 - cdf(Chisq(df), chi2)
293+
p_value = ccdf(Chisq(df), chi2)
294+
295+
min_expected = minimum(expected)
296+
n_low = count(<(5), expected)
297+
warning = min_expected < 5 ?
298+
"$n_low of $k expected cell counts are below 5; the χ² approximation may be unreliable (Cochran's rule of thumb)" :
299+
nothing
197300

198301
return Dict{String,Any}(
199302
"chi_squared" => chi2,
@@ -202,6 +305,8 @@ function chi_square_goodness_of_fit(observed::Vector{Int},
202305
"significant" => p_value < alpha,
203306
"expected" => expected,
204307
"n" => n,
205-
"test_type" => "Chi-square goodness-of-fit"
308+
"test_type" => "Chi-square goodness-of-fit",
309+
"note" => nothing,
310+
"warning" => warning
206311
)
207312
end

src/stats/validation.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,16 @@ function require_probability(v::Real, name::String="value")
8787
return nothing
8888
end
8989

90+
"""
91+
require_nonnegative(x, name::String)
92+
93+
Throw `ArgumentError` unless every element of `x` is `>= 0`.
94+
"""
95+
function require_nonnegative(x, name::String)
96+
all(v -> v >= 0, x) || throw(ArgumentError("$name must not contain negative values"))
97+
return nothing
98+
end
99+
90100
"""
91101
require_square(M::AbstractMatrix, name::String="matrix")
92102

0 commit comments

Comments
 (0)