W1-2: guard partial_correlation against NaN leak on zero-variance input (#40)

## Summary

Fixes the two `must_fix` items from reviewer feedback on W1-2
(`.claude/tasks/prod-readiness/w1-2-stats-degenerate-inputs.md`):
`partial_correlation` still leaked `NaN` on a degenerate `n>=4` input,
and
the degenerate-input test suite never exercised that path.

- **`src/stats/correlation_regression.jl:274` (`partial_correlation`)**
—
  the function guarded `n<4` and the final `denom==0` case, but computed
  `r_xy = cor(x,y)`, `r_xz = cor(x,z)`, `r_yz = cor(y,z)` directly via
  `Statistics.cor()`, which returns `NaN` (not an error) on a
constant/zero-variance vector, and returned that `NaN` straight through
  in the `r_xy`/`r_xz`/`r_yz` fields — violating acceptance criterion #1
  ("No core stat function returns NaN/Inf on degenerate input").
  Replaced the three `cor()` calls with a guarded pairwise-correlation
  helper using the same `den > 0` pattern `pearson_correlation` already
uses, so a constant `x`, `y`, or `z` now yields `r_xy`/`r_xz`/`r_yz` =>
  `nothing` instead of `NaN`, with `r_partial`/`t_stat`/`p_value`
  cascading to `nothing` and a `"note"` explaining why.

- **`test/degenerate_input_test.jl`** — the only existing
  `partial_correlation` degenerate test used `n=3` (the `n<4` early
  return), so the constant-vector/zero-variance case at `n>=4` was never
  tested. Added a new test reproducing the reviewer's exact repro
  (`partial_correlation(fill(2.0,5), [1,2,3,4,5], [2,1,4,3,5])`) that
  asserts `r_xy`, `r_xz`, and `r_partial` are `nothing` and runs the
  `assert_finite_and_serialisable` walk over the full returned `Dict`.

Branch was rebuilt on the latest `main` (which had moved one commit
ahead
with the unrelated P0 guardrail work, #37) via a merge + cherry-pick —
no
force-push, no history rewrite of the already-pushed commit.

## Verification

```
flock /tmp/statistikles-julia.lock -c 'cd <repo> && julia --project=. -e "using Pkg; Pkg.test()"'
```
(WSL Debian login shell.) Full suite green, 0 failures:

| Testset | Pass/Total |
|---|---|
| Statistikles Full Test Suite | 424/424 |
| Statistikles E2E Pipeline Tests | 155/155 |
| Statistikles Property-Based Tests | 3800/3800 |
| Reference Validation (ground truth) | 35/35 |
| Degenerate Input Guards | 943/943 |
| Neural-Boundary Guardrail | 80/80 |

Total: 5437/5437 passed.

## Nothing skipped

Both `must_fix` items from the review are addressed; no other changes
were needed.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

---------

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: hyperpolymath

src/Statistikles.jl

@@ -32,6 +32,7 @@ using Random
 using UUIDs
 
 # --- SYMBOLIC KERNEL: Verified Statistical Methods ---
+include("stats/validation.jl")       # Shared ArgumentError input-validation helpers
 include("stats/descriptive.jl")      # Mean, Median, Variance, Skewness
 include("stats/inferential.jl")      # T-Tests, ANOVA, Chi-Square
 include("stats/correlation_regression.jl")

src/bridge/betlang_bridge.jl

@@ -127,7 +127,8 @@ struct ImpreciseProbability
     lower::Float64
     upper::Float64
     function ImpreciseProbability(lo, hi)
-        @assert 0.0 <= lo <= hi <= 1.0
+        (0.0 <= lo <= hi <= 1.0) || throw(ArgumentError(
+            "ImpreciseProbability bounds must satisfy 0 <= lower <= upper <= 1, got lower=$lo, upper=$hi"))
         new(lo, hi)
     end
 end

src/bridge/typell_levels.jl

@@ -56,7 +56,7 @@ end
 struct Probability
     value::Float64
     function Probability(v::Float64)
-        @assert 0.0 <= v <= 1.0 "Probability must be in [0, 1], got $v"
+        require_probability(v, "Probability")
         new(v)
     end
 end
@@ -126,8 +126,12 @@ struct ModularInt
     modulus::Int
     ModularInt(v::Int, m::Int) = new(mod(v, m), m)
 end
-Base.:+(a::ModularInt, b::ModularInt) = (@assert a.modulus == b.modulus; ModularInt(a.value + b.value, a.modulus))
-Base.:*(a::ModularInt, b::ModularInt) = (@assert a.modulus == b.modulus; ModularInt(a.value * b.value, a.modulus))
+Base.:+(a::ModularInt, b::ModularInt) = a.modulus == b.modulus ?
+    ModularInt(a.value + b.value, a.modulus) :
+    throw(ArgumentError("ModularInt addition requires equal moduli, got $(a.modulus) and $(b.modulus)"))
+Base.:*(a::ModularInt, b::ModularInt) = a.modulus == b.modulus ?
+    ModularInt(a.value * b.value, a.modulus) :
+    throw(ArgumentError("ModularInt multiplication requires equal moduli, got $(a.modulus) and $(b.modulus)"))
 
 # ===========================================================================
 # Level 9-10: Verification provenance types

src/integrations/smtlib_integration.jl

@@ -79,7 +79,7 @@ end
 Verify that p-value correction is monotone: adjusted ≥ raw for all i.
 """
 function smt_verify_correction_monotone(raw_p::Vector{Float64}, adj_p::Vector{Float64})
-    @assert length(raw_p) == length(adj_p)
+    require_equal_length(raw_p, adj_p, "raw_p", "adj_p")
     violations = findall(adj_p .< raw_p .- 1e-15)
     return Dict{String,Any}(
         "monotone" => isempty(violations),

src/stats/assumptions.jl

@@ -3,9 +3,42 @@
 
 function test_normality(data::Vector{Float64})
     n = length(data)
+
+    # DEGENERATE GUARD: the Jarque-Bera statistic needs both skewness
+    # (n >= 3) and kurtosis (n >= 4); below that the moment formulas
+    # divide by zero.
+    if n < 4
+        return Dict{String,Any}(
+            "skewness" => nothing, "kurtosis" => nothing,
+            "jarque_bera" => nothing, "JB_p_value" => nothing,
+            "KS_statistic" => nothing, "n" => n,
+            "normal_skew" => nothing, "normal_kurtosis" => nothing,
+            "overall_assessment" => nothing,
+            "recommendations" => "Insufficient data (need at least 4 observations) to assess normality",
+            "note" => "Skewness/kurtosis-based normality test requires at least 4 observations"
+        )
+    end
+
     sorted = sort(data)
     m = mean(data)
     s = std(data)
+
+    # DEGENERATE GUARD: zero variance makes z-scores 0/0 and the Jarque-Bera
+    # inputs undefined.
+    if s == 0.0
+        return Dict{String,Any}(
+            "skewness" => nothing, "kurtosis" => nothing,
+            "jarque_bera" => nothing, "JB_p_value" => nothing,
+            "KS_statistic" => 0.0, "n" => n,
+            "normal_skew" => nothing, "normal_kurtosis" => nothing,
+            "overall_assessment" => "Data is constant; normality is not meaningfully defined",
+            "recommendations" => n < 30 ?
+                "Small sample — consider non-parametric alternatives" :
+                "Large sample — parametric tests may be robust to mild violations",
+            "note" => "Skewness/kurtosis undefined: zero variance in data"
+        )
+    end
+
     z_scores = (data .- m) ./ s
     skew = (n / ((n - 1) * (n - 2))) * sum(z_scores .^ 3)
     kurt = ((n * (n + 1)) / ((n - 1) * (n - 2) * (n - 3))) * sum(z_scores .^ 4) -
@@ -26,7 +59,8 @@ function test_normality(data::Vector{Float64})
             "Data may violate normality assumption",
         "recommendations" => n < 30 ?
             "Small sample — consider non-parametric alternatives" :
-            "Large sample — parametric tests may be robust to mild violations"
+            "Large sample — parametric tests may be robust to mild violations",
+        "note" => nothing
     )
 end
 

src/stats/correlation_regression.jl

@@ -14,7 +14,20 @@ Uses midranks for tied values. Equivalent to Pearson on ranks.
 """
 function spearman_correlation(x::Vector{Float64}, y::Vector{Float64}; alpha::Float64=0.05)
     n = length(x)
-    @assert length(y) == n "Vectors must be of the same length"
+    require_equal_length(x, y, "x", "y")
+
+    # DEGENERATE GUARD: below n=2, midranks/mean on the ranks either errors
+    # (empty collection) or feeds a negative (n-2) into the significance
+    # formula's sqrt, which throws DomainError. Bail out gracefully instead.
+    if n < 2
+        return Dict{String,Any}(
+            "rho" => nothing, "t_stat" => nothing, "df" => nothing,
+            "p_value" => nothing, "significant" => false,
+            "interpretation" => nothing,
+            "test_type" => "Spearman rank correlation (midranks)",
+            "note" => "Spearman correlation requires at least 2 observations"
+        )
+    end
 
     rx = midranks(x)
     ry = midranks(y)
@@ -51,30 +64,54 @@ LINEAR ASSOCIATION: Computes the Pearson product-moment coefficient.
 """
 function pearson_correlation(x::Vector{Float64}, y::Vector{Float64}; alpha::Float64=0.05)
     n = length(x)
-    @assert length(y) == n "Vectors must be of the same length"
-    
+    require_equal_length(x, y, "x", "y")
+
+    # DEGENERATE GUARD: below n=2 there is nothing to correlate.
+    if n < 2
+        return Dict{String,Any}(
+            "r" => nothing, "r_squared" => nothing, "t_stat" => nothing,
+            "df" => nothing, "p_value" => nothing, "significant" => false,
+            "interpretation" => nothing,
+            "note" => "Pearson correlation requires at least 2 observations"
+        )
+    end
+
     mx, my = mean(x), mean(y)
     num = sum((x .- mx) .* (y .- my))
     den = sqrt(sum((x .- mx).^2) * sum((y .- my).^2))
-    
-    r = num / den
-    
-    # Statistical significance (t-test)
-    # t = r * sqrt((n-2)/(1-r^2))
-    t_stat = r * sqrt((n - 2) / (1 - r^2))
+
+    # DEGENERATE GUARD: den == 0 means zero variance in x or y (r = 0/0).
+    r = den > 0 ? num / den : nothing
+    r_squared = r === nothing ? nothing : r^2
+
+    # Statistical significance (t-test): t = r * sqrt((n-2)/(1-r^2))
+    # DEGENERATE GUARD: needs df > 0 and a non-perfect correlation
+    # (1 - r^2 > 0), else the formula divides by zero.
     df = n - 2
-    p_val = 2 * (1 - cdf(TDist(df), abs(t_stat)))
-    
-    interpretation = abs(r) > 0.7 ? "Strong" : abs(r) > 0.3 ? "Moderate" : "Weak"
-    
+    t_stat = (r !== nothing && df > 0 && (1 - r^2) > 0) ?
+             r * sqrt(df / (1 - r^2)) : nothing
+    p_val = t_stat === nothing ? nothing : 2 * (1 - cdf(TDist(df), abs(t_stat)))
+
+    interpretation = r === nothing ? nothing :
+        abs(r) > 0.7 ? "Strong" : abs(r) > 0.3 ? "Moderate" : "Weak"
+
+    note = if r === nothing
+        "Correlation undefined: zero variance in x or y"
+    elseif t_stat === nothing
+        "Significance test undefined: perfect correlation or insufficient degrees of freedom"
+    else
+        nothing
+    end
+
     return Dict{String,Any}(
         "r" => r,
-        "r_squared" => r^2,
+        "r_squared" => r_squared,
         "t_stat" => t_stat,
         "df" => df,
         "p_value" => p_val,
-        "significant" => p_val < alpha,
-        "interpretation" => interpretation
+        "significant" => p_val === nothing ? false : p_val < alpha,
+        "interpretation" => interpretation,
+        "note" => note
     )
 end
 
@@ -85,24 +122,37 @@ BIVARIATE MODELING: Fits a straight line (y = beta0 + beta1*x).
 """
 function simple_linear_regression(x::Vector{Float64}, y::Vector{Float64}; alpha::Float64=0.05)
     n = length(x)
+    require_equal_length(x, y, "x", "y")
     mx, my = mean(x), mean(y)
-    
-    beta1 = sum((x .- mx) .* (y .- my)) / sum((x .- mx).^2)
+
+    # DEGENERATE GUARD: zero variance in x makes the slope 0/0.
+    ss_xx = sum((x .- mx) .^ 2)
+    if ss_xx == 0
+        return Dict{String,Any}(
+            "intercept" => nothing, "slope" => nothing, "r_squared" => nothing,
+            "n" => n, "test_type" => "Simple Linear Regression",
+            "note" => "Regression undefined: zero variance in x"
+        )
+    end
+
+    beta1 = sum((x .- mx) .* (y .- my)) / ss_xx
     beta0 = my - beta1 * mx
-    
+
     # Predictions and residuals
     y_pred = beta0 .+ beta1 .* x
     residuals = y .- y_pred
     ss_res = sum(residuals.^2)
     ss_tot = sum((y .- my).^2)
-    r2 = 1 - (ss_res / ss_tot)
-    
+    # DEGENERATE GUARD: zero variance in y makes R² a 0/0.
+    r2 = ss_tot > 0 ? 1 - (ss_res / ss_tot) : nothing
+
     return Dict{String,Any}(
         "intercept" => beta0,
         "slope" => beta1,
         "r_squared" => r2,
         "n" => n,
-        "test_type" => "Simple Linear Regression"
+        "test_type" => "Simple Linear Regression",
+        "note" => r2 === nothing ? "R-squared undefined: zero variance in y" : nothing
     )
 end
 
@@ -116,37 +166,57 @@ MULTIVARIATE MODELING: Performs Ordinary Least Squares (OLS) regression.
 """
 function multiple_regression(X::Matrix{Float64}, y::Vector{Float64}; var_names=nothing)
     n, p = size(X)
+    length(y) == n || throw(ArgumentError(
+        "y must have length equal to the number of rows in X (got $(length(y)), expected $n)"))
+
     # Add intercept column
     X_aug = hcat(ones(n), X)
-    
+
     # Normal equation: beta = (XtX) \ Xt * y
     beta = (X_aug' * X_aug) \ (X_aug' * y)
-    
+
     y_pred = X_aug * beta
     residuals = y .- y_pred
     ss_res = sum(residuals.^2)
     ss_tot = sum((y .- mean(y)).^2)
-    
-    r2 = 1 - (ss_res / ss_tot)
-    adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)
-    
-    # Standard errors of coefficients
-    sigma2 = ss_res / (n - p - 1)
-    se = sqrt.(diag(sigma2 * inv(X_aug' * X_aug)))
-    t_stats = beta ./ se
-    p_values = 2 .* (1 .- cdf.(TDist(n - p - 1), abs.(t_stats)))
-    
+
+    # DEGENERATE GUARD: zero variance in y makes R² a 0/0.
+    r2 = ss_tot > 0 ? 1 - (ss_res / ss_tot) : nothing
+    df_resid = n - p - 1
+    adj_r2 = (r2 !== nothing && df_resid > 0) ? 1 - (1 - r2) * (n - 1) / df_resid : nothing
+
     labels = isnothing(var_names) ? ["Intercept"; ["Var$i" for i in 1:p]] : ["Intercept"; var_names]
-    
+
+    # DEGENERATE GUARD: standard errors need positive residual degrees of
+    # freedom (n - p - 1 > 0); otherwise sigma2's denominator is <= 0.
+    note = nothing
+    if r2 !== nothing && df_resid > 0
+        sigma2 = ss_res / df_resid
+        se = sqrt.(diag(sigma2 * inv(X_aug' * X_aug)))
+        t_stats = beta ./ se
+        p_values = 2 .* (1 .- cdf.(TDist(df_resid), abs.(t_stats)))
+        std_errors = Dict(labels[i] => se[i] for i in 1:(p+1))
+        t_stats_dict = Dict(labels[i] => t_stats[i] for i in 1:(p+1))
+        p_values_dict = Dict(labels[i] => p_values[i] for i in 1:(p+1))
+    else
+        std_errors = Dict(labels[i] => nothing for i in 1:(p+1))
+        t_stats_dict = Dict(labels[i] => nothing for i in 1:(p+1))
+        p_values_dict = Dict(labels[i] => nothing for i in 1:(p+1))
+        note = df_resid <= 0 ?
+            "Standard errors undefined: insufficient residual degrees of freedom (n - p - 1 <= 0)" :
+            "Standard errors undefined: zero variance in y"
+    end
+
     return Dict{String,Any}(
         "coefficients" => Dict(labels[i] => beta[i] for i in 1:(p+1)),
-        "std_errors" => Dict(labels[i] => se[i] for i in 1:(p+1)),
-        "t_stats" => Dict(labels[i] => t_stats[i] for i in 1:(p+1)),
-        "p_values" => Dict(labels[i] => p_values[i] for i in 1:(p+1)),
+        "std_errors" => std_errors,
+        "t_stats" => t_stats_dict,
+        "p_values" => p_values_dict,
         "r_squared" => r2,
         "adj_r_squared" => adj_r2,
         "n" => n,
-        "p" => p
+        "p" => p,
+        "note" => note
     )
 end
 
@@ -203,20 +273,61 @@ Removes the effect of confounding variable(s) z.
 """
 function partial_correlation(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64}; alpha::Float64=0.05)
     n = length(x)
-    @assert length(y) == n && length(z) == n
+    require_equal_length(x, y, "x", "y")
+    require_equal_length(x, z, "x", "z")
+
+    # DEGENERATE GUARD: need at least 4 observations (df = n - 3 > 0).
+    if n < 4
+        return Dict{String,Any}(
+            "r_partial" => nothing, "r_xy" => nothing, "r_xz" => nothing, "r_yz" => nothing,
+            "t_stat" => nothing, "df" => nothing, "p_value" => nothing, "significant" => false,
+            "test_type" => "Partial correlation (controlling for z)",
+            "note" => "Partial correlation requires at least 4 observations"
+        )
+    end
+
+    # DEGENERATE GUARD: Statistics.cor() returns NaN (not an error) for a
+    # constant/zero-variance input, which would otherwise leak straight into
+    # r_xy/r_xz/r_yz. Compute each pairwise correlation the same guarded way
+    # pearson_correlation does (den == 0 => nothing) so no NaN can enter the
+    # partial-correlation formula below.
+    safe_cor(a::Vector{Float64}, b::Vector{Float64}) = begin
+        ma, mb = mean(a), mean(b)
+        den = sqrt(sum((a .- ma) .^ 2) * sum((b .- mb) .^ 2))
+        den > 0 ? sum((a .- ma) .* (b .- mb)) / den : nothing
+    end
 
-    r_xy = cor(x, y)
-    r_xz = cor(x, z)
-    r_yz = cor(y, z)
+    r_xy = safe_cor(x, y)
+    r_xz = safe_cor(x, z)
+    r_yz = safe_cor(y, z)
 
     # Partial r formula: r_xy.z = (r_xy - r_xz * r_yz) / sqrt((1 - r_xz²)(1 - r_yz²))
-    denom = sqrt((1 - r_xz^2) * (1 - r_yz^2))
-    r_partial = denom > 0 ? (r_xy - r_xz * r_yz) / denom : 0.0
+    # DEGENERATE GUARD: undefined if any pairwise correlation above is
+    # undefined (zero variance in x, y, or z), or if the residual variance
+    # denominator is zero.
+    r_partial = if r_xy === nothing || r_xz === nothing || r_yz === nothing
+        nothing
+    else
+        denom = sqrt((1 - r_xz^2) * (1 - r_yz^2))
+        denom > 0 ? (r_xy - r_xz * r_yz) / denom : nothing
+    end
 
-    # Significance test
+    # Significance test. DEGENERATE GUARD: needs df > 0 and a non-perfect
+    # partial correlation (1 - r_partial^2 > 0).
     df = n - 3
-    t_stat = r_partial * sqrt(df / (1 - r_partial^2))
-    p_val = df > 0 ? 2 * (1 - cdf(TDist(df), abs(t_stat))) : 1.0
+    t_stat = (r_partial !== nothing && df > 0 && (1 - r_partial^2) > 0) ?
+             r_partial * sqrt(df / (1 - r_partial^2)) : nothing
+    p_val = t_stat === nothing ? nothing : 2 * (1 - cdf(TDist(df), abs(t_stat)))
+
+    note = if r_xy === nothing || r_xz === nothing || r_yz === nothing
+        "Partial correlation undefined: zero variance in x, y, or z (constant input)"
+    elseif r_partial === nothing
+        "Partial correlation undefined: zero residual variance after controlling for z"
+    elseif t_stat === nothing
+        "Significance test undefined: perfect partial correlation or insufficient degrees of freedom"
+    else
+        nothing
+    end
 
     return Dict{String,Any}(
         "r_partial" => r_partial,
@@ -226,8 +337,9 @@ function partial_correlation(x::Vector{Float64}, y::Vector{Float64}, z::Vector{F
         "t_stat" => t_stat,
         "df" => df,
         "p_value" => p_val,
-        "significant" => p_val < alpha,
-        "test_type" => "Partial correlation (controlling for z)"
+        "significant" => p_val === nothing ? false : p_val < alpha,
+        "test_type" => "Partial correlation (controlling for z)",
+        "note" => note
     )
 end