Statistical tests

use statistics.qweruse stat_ttest, stat_anova from statistics      // Quillo 0.6: only these

statrs = "0.17"

// A clinic's and a shop's questions answered with tests: did the new treatment lower blood// pressure (paired t-test)? do three page designs sell differently (ANOVA)? does the choice of// plan depend on the country (chi-square)? and what drives a flat's price (multiple regression)?////   quillo run statistics.qwe      then open http://localhost:8080/tests and http://localhost:8080/regressionuse "statistics.qwer"let BEFORE = [148.0, 152.0, 139.0, 160.0, 155.0, 147.0, 151.0, 158.0]let AFTER = [141.0, 147.0, 138.0, 150.0, 149.0, 145.0, 144.0, 151.0]route GET /tests = {    paired: stat_ttest("paired", BEFORE, AFTER, 0.0, "two")?,    welch: stat_ttest("welch", [12.1, 14.3, 11.8, 13.5, 12.9], [15.2, 16.8, 14.9, 17.1], 0.0, "less")?,    anova: stat_anova([[23.0, 25.0, 21.0, 24.0], [30.0, 28.0, 33.0, 29.0], [24.0, 26.0, 22.0, 25.0]])?,    plans: stat_chi2_table([[120.0, 90.0, 40.0], [80.0, 110.0, 60.0]], false)?,    dice: stat_chi2_fit([8.0, 9.0, 12.0, 11.0, 6.0, 14.0], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0])?,    waits: stat_mann_whitney([3.1, 4.5, 2.8, 5.2, 3.9], [6.1, 7.4, 5.9, 8.2, 6.6, 7.0])?,    ads: stat_correlation([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [2.1, 3.9, 6.2, 7.8, 10.1, 12.2], false)?,    ci: stat_mean_ci(AFTER, 0.95)?,}let FLATS = [    { price: 210.0, m2: 52.0, rooms: 2.0, floor: 1.0 }, { price: 265.0, m2: 68.0, rooms: 3.0, floor: 2.0 }, { price: 180.0, m2: 45.0, rooms: 1.0, floor: 0.0 },    { price: 330.0, m2: 85.0, rooms: 3.0, floor: 5.0 }, { price: 295.0, m2: 74.0, rooms: 3.0, floor: 3.0 }, { price: 240.0, m2: 60.0, rooms: 2.0, floor: 4.0 },    { price: 390.0, m2: 101.0, rooms: 4.0, floor: 2.0 }, { price: 205.0, m2: 50.0, rooms: 2.0, floor: 0.0 }, { price: 350.0, m2: 92.0, rooms: 4.0, floor: 1.0 },]route GET /regression = stat_regression(FLATS, "price", ["m2", "rooms", "floor"])?

//! Statistical tests and regression beyond Quillo's `stats`: one-sample, paired, Student's and Welch's t-tests, one-way ANOVA, chi-square tests of independence and goodness of fit, Mann-Whitney U, correlation with its p-value, confidence intervals, and multiple linear regression with standard errors, t and p values, R² and predictions.//! version: 1.0.0//!//! [dependencies]//! statrs = "0.17"use serde_json::{json, Value};use statrs::distribution::{ChiSquared, ContinuousCDF, FisherSnedecor, Normal, StudentsT};fn r(x: f64) -> Value {    if x.is_finite() { json!((x * 1e6).round() / 1e6) } else { Value::Null }}fn mean(x: &[f64]) -> f64 {    x.iter().sum::<f64>() / x.len() as f64}fn var(x: &[f64]) -> f64 {    let m = mean(x);    x.iter().map(|v| (v - m).powi(2)).sum::<f64>() / (x.len() as f64 - 1.0)}fn t_p(t: f64, df: f64, tail: &str) -> Result<f64, String> {    let d = StudentsT::new(0.0, 1.0, df).map_err(|e| e.to_string())?;    Ok(match tail {        "less" => d.cdf(t),        "greater" => 1.0 - d.cdf(t),        _ => 2.0 * (1.0 - d.cdf(t.abs())),    })}fn t_crit(df: f64, conf: f64) -> Result<f64, String> {    Ok(StudentsT::new(0.0, 1.0, df).map_err(|e| e.to_string())?.inverse_cdf(1.0 - (1.0 - conf) / 2.0))}fn check(x: &[f64], name: &str, min: usize) -> Result<(), String> {    if x.len() < min {        return Err(format!("{} needs at least {} values", name, min));    }    if x.iter().any(|v| !v.is_finite()) {        return Err(format!("{} has a value that is not a number", name));    }    Ok(())}/// A t-test. `kind`: "one" (is `a`'s mean `mu`?), "paired" (`a` and `b` measured on the same subjects), "welch" (two groups, variances may differ; the safe default) or "student" (equal variances)./// `tail`: "two", "less" or "greater". Gives `{ t, df, p, mean_difference, ci (95% of the difference), cohen_d, significant (p < 0.05) }`.pub fn stat_ttest(kind: &str, a: Vec<f64>, b: Vec<f64>, mu: f64, tail: &str) -> Result<serde_json::Value, String> {    check(&a, "a", 2)?;    let (t, df, diff, se, d) = match kind {        "one" => {            let se = (var(&a) / a.len() as f64).sqrt();            let diff = mean(&a) - mu;            (diff / se, a.len() as f64 - 1.0, diff, se, diff / var(&a).sqrt())        }        "paired" => {            check(&b, "b", 2)?;            if a.len() != b.len() {                return Err("paired samples have the same length".into());            }            let dd: Vec<f64> = a.iter().zip(&b).map(|(x, y)| x - y).collect();            let se = (var(&dd) / dd.len() as f64).sqrt();            let diff = mean(&dd) - mu;            (diff / se, dd.len() as f64 - 1.0, diff, se, diff / var(&dd).sqrt())        }        "welch" | "student" => {            check(&b, "b", 2)?;            let (na, nb) = (a.len() as f64, b.len() as f64);            let (va, vb) = (var(&a), var(&b));            let diff = mean(&a) - mean(&b) - mu;            let pooled = (((na - 1.0) * va + (nb - 1.0) * vb) / (na + nb - 2.0)).sqrt();            if kind == "student" {                let se = pooled * (1.0 / na + 1.0 / nb).sqrt();                (diff / se, na + nb - 2.0, diff, se, diff / pooled)            } else {                let se = (va / na + vb / nb).sqrt();                let df = (va / na + vb / nb).powi(2) / ((va / na).powi(2) / (na - 1.0) + (vb / nb).powi(2) / (nb - 1.0));                (diff / se, df, diff, se, diff / pooled)            }        }        other => return Err(format!("kind {:?}: one, paired, welch or student", other)),    };    if !se.is_finite() || se == 0.0 {        return Err("the values do not vary: a t-test needs some spread".into());    }    let p = t_p(t, df, tail)?;    let tc = t_crit(df, 0.95)?;    Ok(json!({ "t": r(t), "df": r(df), "p": r(p), "mean_difference": r(diff), "ci": [r(diff - tc * se), r(diff + tc * se)], "cohen_d": r(d), "significant": p < 0.05 }))}/// One-way ANOVA: do these groups' means differ? Gives `{ f, df_between, df_within, p, eta_squared, significant, means }`.pub fn stat_anova(groups: Vec<Vec<f64>>) -> Result<serde_json::Value, String> {    if groups.len() < 2 {        return Err("ANOVA compares two or more groups".into());    }    for (i, g) in groups.iter().enumerate() {        check(g, &format!("group {}", i + 1), 2)?;    }    let all: Vec<f64> = groups.iter().flatten().copied().collect();    let gm = mean(&all);    let ssb: f64 = groups.iter().map(|g| g.len() as f64 * (mean(g) - gm).powi(2)).sum();    let ssw: f64 = groups.iter().map(|g| {        let m = mean(g);        g.iter().map(|v| (v - m).powi(2)).sum::<f64>()    }).sum();    let (dfb, dfw) = (groups.len() as f64 - 1.0, all.len() as f64 - groups.len() as f64);    let f = (ssb / dfb) / (ssw / dfw);    let p = 1.0 - FisherSnedecor::new(dfb, dfw).map_err(|e| e.to_string())?.cdf(f);    Ok(json!({ "f": r(f), "df_between": dfb, "df_within": dfw, "p": r(p), "eta_squared": r(ssb / (ssb + ssw)), "significant": p < 0.05, "means": groups.iter().map(|g| r(mean(g))).collect::<Vec<_>>() }))}/// Chi-square test of independence on a table of counts (rows by columns: did choice depend on group?). Gives `{ chi2, df, p, cramers_v, expected, significant }`; with `yates` for a 2 × 2 table.pub fn stat_chi2_table(table: Vec<Vec<f64>>, yates: bool) -> Result<serde_json::Value, String> {    let rows = table.len();    let cols = table.first().map_or(0, |r| r.len());    if rows < 2 || cols < 2 || table.iter().any(|r| r.len() != cols) {        return Err("a table of at least 2 rows and 2 columns, all rows as long".into());    }    if table.iter().flatten().any(|v| *v < 0.0 || !v.is_finite()) {        return Err("counts are 0 or more".into());    }    let total: f64 = table.iter().flatten().sum();    let rs: Vec<f64> = table.iter().map(|r| r.iter().sum()).collect();    let cs: Vec<f64> = (0..cols).map(|j| table.iter().map(|r| r[j]).sum()).collect();    let mut chi = 0.0;    let mut expected = vec![vec![0.0; cols]; rows];    for i in 0..rows {        for j in 0..cols {            let e = rs[i] * cs[j] / total;            expected[i][j] = e;            if e > 0.0 {                let d = (table[i][j] - e).abs() - if yates && rows == 2 && cols == 2 { 0.5 } else { 0.0 };                chi += d.max(0.0).powi(2) / e;            }        }    }    let df = ((rows - 1) * (cols - 1)) as f64;    let p = 1.0 - ChiSquared::new(df).map_err(|e| e.to_string())?.cdf(chi);    let v = (chi / (total * ((rows.min(cols) - 1) as f64))).sqrt();    let low = expected.iter().flatten().filter(|e| **e < 5.0).count();    Ok(json!({ "chi2": r(chi), "df": df, "p": r(p), "cramers_v": r(v), "significant": p < 0.05,        "expected": expected.iter().map(|row| row.iter().map(|e| r(*e)).collect::<Vec<_>>()).collect::<Vec<_>>(),        "warning": if low > 0 { json!(format!("{} expected counts are under 5: the test is rough", low)) } else { Value::Null } }))}/// Chi-square goodness of fit: do `observed` counts fit `expected` proportions (or counts; scaled to the same total)? Gives `{ chi2, df, p, significant }`.pub fn stat_chi2_fit(observed: Vec<f64>, expected: Vec<f64>) -> Result<serde_json::Value, String> {    if observed.len() < 2 || observed.len() != expected.len() {        return Err("observed and expected, as long as each other, at least 2".into());    }    let (to, te): (f64, f64) = (observed.iter().sum(), expected.iter().sum());    if te <= 0.0 {        return Err("expected values add up to more than 0".into());    }    let chi: f64 = observed.iter().zip(&expected).map(|(o, e)| {        let e = e / te * to;        if e > 0.0 { (o - e).powi(2) / e } else { 0.0 }    }).sum();    let df = observed.len() as f64 - 1.0;    let p = 1.0 - ChiSquared::new(df).map_err(|e| e.to_string())?.cdf(chi);    Ok(json!({ "chi2": r(chi), "df": df, "p": r(p), "significant": p < 0.05 }))}/// Ranks with ties sharing their average rank.fn ranks(v: &[f64]) -> Vec<f64> {    let mut idx: Vec<usize> = (0..v.len()).collect();    idx.sort_by(|a, b| v[*a].partial_cmp(&v[*b]).unwrap());    let mut out = vec![0.0; v.len()];    let mut i = 0;    while i < idx.len() {        let mut j = i;        while j + 1 < idx.len() && v[idx[j + 1]] == v[idx[i]] {            j += 1;        }        let rank = (i + j) as f64 / 2.0 + 1.0;        for k in i..=j {            out[idx[k]] = rank;        }        i = j + 1;    }    out}/// Mann-Whitney U (Wilcoxon rank-sum): do two groups differ when the data are not normal (times, ratings)? Normal approximation with ties corrected. Gives `{ u, z, p, significant, effect (rank-biserial) }`.pub fn stat_mann_whitney(a: Vec<f64>, b: Vec<f64>) -> Result<serde_json::Value, String> {    check(&a, "a", 1)?;    check(&b, "b", 1)?;    let (n1, n2) = (a.len() as f64, b.len() as f64);    let all: Vec<f64> = a.iter().chain(&b).copied().collect();    let rk = ranks(&all);    let r1: f64 = rk[..a.len()].iter().sum();    let u1 = r1 - n1 * (n1 + 1.0) / 2.0;    let u = u1.min(n1 * n2 - u1);    let n = n1 + n2;    let mut ties = 0.0;    let mut sorted = all.clone();    sorted.sort_by(|x, y| x.partial_cmp(y).unwrap());    let mut i = 0;    while i < sorted.len() {        let mut j = i;        while j + 1 < sorted.len() && sorted[j + 1] == sorted[i] {            j += 1;        }        let t = (j - i + 1) as f64;        ties += t * t * t - t;        i = j + 1;    }    let sd = (n1 * n2 / 12.0 * ((n + 1.0) - ties / (n * (n - 1.0)))).sqrt();    let z = if sd > 0.0 { (u1 - n1 * n2 / 2.0) / sd } else { 0.0 };    let p = 2.0 * (1.0 - Normal::new(0.0, 1.0).unwrap().cdf(z.abs()));    Ok(json!({ "u": r(u), "u1": r(u1), "z": r(z), "p": r(p), "significant": p < 0.05, "effect": r(1.0 - 2.0 * u / (n1 * n2)) }))}/// Pearson's (or, with `ranked`, Spearman's) correlation with its p-value and 95% interval: `{ r, p, ci, n, significant }`.pub fn stat_correlation(x: Vec<f64>, y: Vec<f64>, ranked: bool) -> Result<serde_json::Value, String> {    check(&x, "x", 3)?;    if x.len() != y.len() {        return Err("x and y have the same length".into());    }    let (x, y) = if ranked { (ranks(&x), ranks(&y)) } else { (x, y) };    let (mx, my) = (mean(&x), mean(&y));    let sxy: f64 = x.iter().zip(&y).map(|(a, b)| (a - mx) * (b - my)).sum();    let sxx: f64 = x.iter().map(|a| (a - mx).powi(2)).sum();    let syy: f64 = y.iter().map(|b| (b - my).powi(2)).sum();    if sxx == 0.0 || syy == 0.0 {        return Err("a list that does not vary has no correlation".into());    }    let rr = sxy / (sxx * syy).sqrt();    let n = x.len() as f64;    let t = rr * ((n - 2.0) / (1.0 - rr * rr).max(1e-300)).sqrt();    let p = t_p(t, n - 2.0, "two")?;    // Fisher's z for the interval.    let ci = if n > 3.0 && rr.abs() < 1.0 {        let z = rr.atanh();        let se = 1.0 / (n - 3.0).sqrt();        json!([r((z - 1.96 * se).tanh()), r((z + 1.96 * se).tanh())])    } else {        Value::Null    };    Ok(json!({ "r": r(rr), "p": r(p), "ci": ci, "n": n as i64, "significant": p < 0.05 }))}/// The confidence interval of a mean (t-based): `{ mean, low, high, margin, n }` at `confidence` (0.95).pub fn stat_mean_ci(x: Vec<f64>, confidence: f64) -> Result<serde_json::Value, String> {    check(&x, "x", 2)?;    let conf = if confidence > 0.0 && confidence < 1.0 { confidence } else { 0.95 };    let m = mean(&x);    let margin = t_crit(x.len() as f64 - 1.0, conf)? * (var(&x) / x.len() as f64).sqrt();    Ok(json!({ "mean": r(m), "low": r(m - margin), "high": r(m + margin), "margin": r(margin), "n": x.len() }))}/// Solves A·x = b by Gauss-Jordan with pivoting; also gives A's inverse.fn solve(mut a: Vec<Vec<f64>>) -> Option<Vec<Vec<f64>>> {    let n = a.len();    let mut inv: Vec<Vec<f64>> = (0..n).map(|i| (0..n).map(|j| if i == j { 1.0 } else { 0.0 }).collect()).collect();    for c in 0..n {        let p = (c..n).max_by(|x, y| a[*x][c].abs().partial_cmp(&a[*y][c].abs()).unwrap())?;        if a[p][c].abs() < 1e-12 {            return None;        }        a.swap(c, p);        inv.swap(c, p);        let d = a[c][c];        for j in 0..n {            a[c][j] /= d;            inv[c][j] /= d;        }        for i in 0..n {            if i != c {                let f = a[i][c];                for j in 0..n {                    a[i][j] -= f * a[c][j];                    inv[i][j] -= f * inv[c][j];                }            }        }    }    Some(inv)}/// Multiple linear regression by least squares: `rows` are maps (or lists), `y` the field predicted, `x` the predictors' fields (or column numbers)./// Gives `{ coefficients: [{ name, estimate, std_error, t, p }], r2, adj_r2, f, f_p, residual_se, n }`, the intercept first.pub fn stat_regression(rows: Vec<serde_json::Value>, y: serde_json::Value, x: Vec<serde_json::Value>) -> Result<serde_json::Value, String> {    let get = |row: &Value, k: &Value| -> Option<f64> {        match k {            Value::String(s) => row.get(s.as_str()).and_then(Value::as_f64),            Value::Number(n) => row.get(n.as_u64()? as usize).and_then(Value::as_f64),            _ => None,        }    };    let p = x.len() + 1;    let mut xs: Vec<Vec<f64>> = vec![];    let mut ys: Vec<f64> = vec![];    for (i, row) in rows.iter().enumerate() {        let yv = get(row, &y).ok_or_else(|| format!("row {}: no number for {}", i + 1, y))?;        let mut xr = vec![1.0];        for k in &x {            xr.push(get(row, k).ok_or_else(|| format!("row {}: no number for {}", i + 1, k))?);        }        xs.push(xr);        ys.push(yv);    }    let n = xs.len();    if n <= p {        return Err(format!("{} predictors need more than {} rows", x.len(), p));    }    let xtx: Vec<Vec<f64>> = (0..p).map(|i| (0..p).map(|j| xs.iter().map(|r| r[i] * r[j]).sum()).collect()).collect();    let xty: Vec<f64> = (0..p).map(|i| xs.iter().zip(&ys).map(|(r, yv)| r[i] * yv).sum()).collect();    let inv = solve(xtx).ok_or("the predictors are collinear (one is a mix of the others)")?;    let beta: Vec<f64> = (0..p).map(|i| (0..p).map(|j| inv[i][j] * xty[j]).sum()).collect();    let fitted: Vec<f64> = xs.iter().map(|r| r.iter().zip(&beta).map(|(a, b)| a * b).sum()).collect();    let my = mean(&ys);    let sse: f64 = ys.iter().zip(&fitted).map(|(a, b)| (a - b).powi(2)).sum();    let sst: f64 = ys.iter().map(|a| (a - my).powi(2)).sum();    let df = (n - p) as f64;    let s2 = sse / df;    let r2 = if sst > 0.0 { 1.0 - sse / sst } else { 1.0 };    let names: Vec<String> = std::iter::once("intercept".to_string()).chain(x.iter().map(|k| match k { Value::String(s) => s.clone(), o => o.to_string() })).collect();    let mut coefs = vec![];    for i in 0..p {        let se = (s2 * inv[i][i]).sqrt();        let t = beta[i] / se;        coefs.push(json!({ "name": names[i], "estimate": r(beta[i]), "std_error": r(se), "t": r(t), "p": r(t_p(t, df, "two")?) }));    }    let k = (p - 1) as f64;    let f = if k > 0.0 && sse > 0.0 { ((sst - sse) / k) / s2 } else { f64::NAN };    let f_p = if f.is_finite() { 1.0 - FisherSnedecor::new(k, df).map_err(|e| e.to_string())?.cdf(f) } else { f64::NAN };    Ok(json!({ "coefficients": coefs, "r2": r(r2), "adj_r2": r(1.0 - (1.0 - r2) * (n as f64 - 1.0) / df), "f": r(f), "f_p": r(f_p), "residual_se": r(s2.sqrt()), "n": n }))}