`mod stats`¶

module stats¶

Statistical functions: chi-squared test with Yates correction, Bonferroni multiple testing correction, and group bias.

Functions

fn bayes_factor_2x2(n_g1: u32, n_g2: u32, total_g1: u32, total_g2: u32) -> f64¶: Bayes Factor for association in a 2x2 contingency table. BF > 1: evidence for association. BF > 10: strong evidence. Uses the version-1 uniform Beta prior compatibility profile.

fn bayes_factor_2x2_with_model(n_g1: u32, n_g2: u32, total_g1: u32, total_g2: u32, model: &BayesFactorModel) -> Result<f64, BayesianModelError>¶: Bayes Factor with independently configurable Beta priors.

fn benjamini_hochberg(p_values: &[f64]) -> Vec<f64>¶: Apply Benjamini-Hochberg FDR correction to a vector of p-values. Returns adjusted p-values (q-values). Controls FDR at the given level.

fn benjamini_hochberg_weighted(p_values: &[(f64, u64)]) -> Vec<f64>¶

Apply Benjamini-Hochberg FDR correction to compressed p-value groups.

Each tuple contains a p-value and the number of marker tests represented by that value. Entries with zero weight do not belong to the testing universe and receive an adjusted p-value of one.

fn bonferroni_correct(p: f64, n_markers: u64) -> f64¶: Apply Bonferroni correction to a p-value.

fn chi_squared_p(chi_sq: f64) -> f64¶

P-value for a chi-squared statistic with df=1.

Uses the exact identity: for df=1, the chi-squared CDF is: P(chi2) = erf(sqrt(chi2/2))
so the p-value is:: p = 1 - P(chi2) = erfc(sqrt(chi2/2))

This replaces the full regularized gamma function with a single libm erfc call. Derived via SymPy:

gamma(1/2, x) = sqrt(pi) * erf(sqrt(x))
Gamma(1/2) = sqrt(pi)
P(1/2, x) = erf(sqrt(x))
p = erfc(sqrt(chi2/2))

fn chi_squared_yates(n_group1: u32, n_group2: u32, total_group1: u32, total_group2: u32) -> f64¶

Chi-squared statistic with Yates continuity correction for a 2x2 table.

Implements the shortcut formula:

chi2 = N * (|ad - bc| - N/2)^2 / (a+b)(c+d)(a+c)(b+d)

where the contingency table is:

====== ============== =============
\      marker present marker absent
====== ============== =============
group1 n_group1       total1 - n1
group2 n_group2       total2 - n2
====== ============== =============

fn empirical_bayes_em(group_counts: &[(u32, u32)], total_g1: u32, total_g2: u32, p_sex: f64, max_iter: usize) -> (f64, Vec<f64>)¶: Empirical Bayes EM: estimate pi (fraction of sex-linked markers) from data. Returns (pi, posteriors) after convergence. group_counts: Vec of (n_g1, n_g2) for each marker.

fn erfc_panelled(t: f64) -> f64¶

erfc(t) for t >= 0 from a panelled minimax table, without libm.

[0, 6) is cut into 24 panels of width 1/4 and each panel carries a degree-14 fit in the centred variable u = t - c, so |u| <= 1/8 and every monomial the Horner loop touches stays O(1). That is what keeps the evaluation error near the approximation error; a single wide fit in the raw variable does not, because its intermediates run to 6^40.

fpminimax bounds the fitted polynomials at 1.1e-16 relative for t < 4.75; evaluating a degree-14 Horner adds a few ulp, and the measured figure against libm is 4.4e-16. Past t = 4.75 it grows to 2.8e-15 by t = 6, where erfc is already under the 1e-16 floor rsx reports p-values at, and past t = 6 the table returns 0 because what it drops is smaller than that floor. scripts/sollya/erfc_panels.sollya generates the table and prints the fitted bound per panel; tests/test_precision.rs measures the evaluated error, which is the number that matters and is never the one the generator reports.

Provided for callers that cannot rely on the platform erfc. chi_squared_p uses libm, which is faster on the hosts we benchmark.

fn find_median(data: &mut [u16]) -> f64¶

Find the mathematical median of a mutable slice (partially reorders in-place).

Uses order-statistic selection (select_nth_unstable) at index len/2 instead of a full sort. Odd lengths return that middle value. Even lengths average it with the maximum of the lower partition, which is the adjacent order statistic. Selection is O(n) average vs O(n log n) sort; the rank identities are validated in proofs/lean/MedianSelect/MedianSelect.lean and scripts/sympy/median_select_proof.py.

fn fisher_exact(n_g1: u32, n_g2: u32, total_g1: u32, total_g2: u32) -> f64¶

Fisher’s exact test for a 2×2 table (two-sided, probability method).

Sums hypergeometric probabilities of all tables with probability less than or equal to the observed table (within a small log-space tolerance). This matches the common two-sided definition used by R fisher.test / SciPy fisher_exact(..., alternative="two-sided") density method — not a one-sided greater/less tail.

Returns 1.0 (clamped floor 1e-16 applied only for tiny tails) when the input is not a valid contingency table (n_g* > total_g*).

fn g_test(n_g1: u32, n_g2: u32, total_g1: u32, total_g2: u32) -> f64¶: G-test (log-likelihood ratio) for 2x2 table. Better asymptotic properties than chi-squared. Returns 1.0 if present counts exceed group totals (invalid table).

fn group_bias(n_group1: u32, total_group1: u32, n_group2: u32, total_group2: u32) -> f64¶: Group bias: difference in marker frequency between two groups. Ranges from -1.0 (only in group2) to +1.0 (only in group1). Returns 0.0 if either group has zero individuals (undefined frequency).

fn logistic_regression(x: &[f64], y: &[f64], n: usize, p: usize, max_iter: usize) -> Vec<f64>¶: Logistic regression: fit y ~ X using IRLS (Newton-Raphson). X: n x p design matrix (row-major), y: n binary outcomes (0/1). Returns coefficient vector beta (length p).

fn p_association(n_group1: u32, n_group2: u32, total_group1: u32, total_group2: u32) -> f64¶: Compute p-value of association with group using chi-squared test with Yates correction. Matches C++ get_p_association exactly.

fn posterior_sex_linked(n_g1: u32, n_g2: u32, total_g1: u32, total_g2: u32, pi: f64, p_sex: f64) -> f64¶

Posterior probability under the compatibility parameter surface.

Calculation paths that expose the complete model use posterior_sex_linked_with_model so the null prevalence and directional mixture weight are explicit.

fn posterior_sex_linked_with_model(n_g1: u32, n_g2: u32, total_g1: u32, total_g2: u32, model: &DirectionalModel) -> f64¶: Posterior under a directional mixture whose calculation inputs are explicit.