The statistics layer

use malevich::stat::Bins;
use malevich::{Bars, Plot, Scale};
use super::penguins;
// The histogram, spelled out: Bins::auto chooses the geometry, Bars::spans draws the counts.
let mass = penguins().mass;
let mut bins = Bins::auto(&mass, 60).expect("finite data");
for value in &mass {
    bins.add(*value);
}
let counts: Vec<f64> = bins.counts().iter().map(|&c| c as f64).collect();
Plot::new()
    .layer(Bars::spans(bins.start(), bins.width(), counts))
    .y_scale(Scale::Integer)
    .title("Bins::auto + Bars::spans + Scale::Integer == hist")
    .x_label("body mass, g")

use malevich::stat::Normalization;
use malevich::{hist_with, HistogramOptions};
let mut unit = super::noise(3);
// Latencies: log-normal, the long tail every service has.
let latency: Vec<f64> = (0..4000).map(|_| (0.4 * super::gaussian(&mut unit) + 3.2).exp()).collect();
hist_with(latency, HistogramOptions::new(40).normalization(Normalization::Percent).cumulative(true))
    .expect("valid options")
    .title("cumulative share of requests")
    .x_label("ms")
    .y_label("%")

use malevich::stat::{calendar_bins, TimeUnit};
use malevich::{Bars, Plot};
use super::day_stamp;
let mut unit = super::noise(29);
// Incidents over a year: counted per calendar month, each bar drawn between
// its month's true edges on the time axis (February is shorter).
let stamps: Vec<f64> = (0..300).map(|_| day_stamp(2025, 1, 1) + unit() * 365.0 * 86_400.0 * (0.4 + 0.6 * unit())).collect();
let months = calendar_bins(&stamps, TimeUnit::Month).expect("data");
let counts: Vec<f64> = months.counts().iter().map(|&c| c as f64).collect();
Plot::new()
    .layer(Bars::intervals(months.starts().to_vec(), months.ends().to_vec(), counts))
    .time_x()
    .title("stat::calendar_bins — incidents per month")

use malevich::stat::{kde_with, Bandwidth, KdeOptions};
use malevich::{Line, Plot};
use super::penguins;
// The bandwidth is the one honest knob a density has: the same 342 masses
// at Silverman's rule, at half of it, and at twice it.
let mass = penguins().mass;
let mut plot = Plot::new();
for (scale, label) in [(0.5, "½ Silverman"), (1.0, "Silverman (default)"), (2.0, "2× Silverman")] {
    let options = KdeOptions::new().bandwidth(Bandwidth::Scale(scale));
    let (x, y) = kde_with(&mass, 200, options).expect("valid").expect("data");
    plot = plot.layer(Line::xy(x, y).label(label));
}
plot.title("kernel density, three bandwidths").x_label("body mass, g")

use malevich::stat::Whiskers;
use malevich::{box_plot_with, BoxOptions};
use super::penguins;
let p = penguins();
let groups = ["Adelie", "Chinstrap", "Gentoo"].map(|s| p.of(s, &p.mass));
// Whiskers to the 5th and 95th percentiles instead of Tukey's 1.5 IQR.
box_plot_with(["Adelie", "Chinstrap", "Gentoo"], groups, BoxOptions::new().whiskers(Whiskers::Percentiles(0.05, 0.95)))
    .expect("valid")
    .title("body mass, whiskers at p5 and p95")
    .y_label("g")

use super::penguins;
let p = penguins();
let groups = ["Adelie", "Chinstrap", "Gentoo"].map(|s| p.of(s, &p.flipper));
// Moments, quantiles, and extremes as a stat table: text on two band axes.
malevich::describe(["Adelie", "Chinstrap", "Gentoo"], groups).title("flipper length, mm")

use malevich::stat::{Reducer, Window};
use malevich::{Line, Plot};
let mut unit = super::noise(5);
// Latency samples with a spiky tail; a rolling mean and a rolling p95 over 50.
let samples: Vec<f64> = (0..400).map(|i| {
    let base = 40.0 + 10.0 * (f64::from(i) * 0.03).sin();
    if unit() < 0.06 { base + 60.0 * unit() } else { base + 6.0 * super::gaussian(&mut unit) }
}).collect();
let window = Window::new(50);
Plot::new()
    .layer(Line::y(samples.clone()).label("samples"))
    .layer(Line::y(window.mean(&samples)).label("rolling mean"))
    .layer(Line::y(window.reduce(&samples, Reducer::Percentile(0.95))).label("rolling p95"))
    .title("Window: one reducer vocabulary")
    .y_label("ms")

use malevich::stat::{Window, WindowAnchor};
use malevich::{Line, Plot};
let mut unit = super::noise(8);
let signal: Vec<f64> = (0..200).map(|i| if (60..70).contains(&i) { 10.0 } else { 2.0 + super::gaussian(&mut unit) * 0.4 }).collect();
// A trailing window lags the step; a centered one straddles it; `strict`
// gaps the positions whose window is incomplete instead of guessing.
Plot::new()
    .layer(Line::y(signal.clone()).label("signal"))
    .layer(Line::y(Window::new(21).mean(&signal)).label("trailing"))
    .layer(Line::y(Window::new(21).anchor(WindowAnchor::Middle).strict().mean(&signal)).label("centered, strict"))
    .title("window anchors")

use malevich::stat::{binned, Bins, Reducer};
use malevich::{Bars, Line, Plot, Points};
let mut unit = super::noise(19);
// Confidence versus accuracy: 0/1 outcomes, binned by confidence, mean per bin.
let confidence: Vec<f64> = (0..2000).map(|_| unit()).collect();
let correct: Vec<f64> = confidence.iter().map(|c| if unit() < c * 0.8 + 0.1 { 1.0 } else { 0.0 }).collect();
let bins = Bins::new(0.0, 0.1, 10);
let accuracy = binned(&confidence, &correct, &bins, Reducer::Mean);
let centers: Vec<f64> = (0..10).map(|i| 0.05 + 0.1 * i as f64).collect();
Plot::new()
    .layer(Bars::spans(0.0, 0.1, accuracy.clone()).label("accuracy per bin"))
    .layer(Points::xy(centers, accuracy))
    .layer(Line::xy(vec![0.0, 1.0], vec![0.0, 1.0]).label("perfect calibration"))
    .title("stat::binned with Reducer::Mean")
    .x_label("confidence")
    .y_label("accuracy")

use malevich::stat::{Agg, Reducer};
use malevich::{Bars, Plot};
use super::penguins;
// Group by a key, reduce each group: the median mass per species.
let p = penguins();
let (species, medians) = Agg::by(p.species, p.mass).reduce(Reducer::Median);
Plot::new()
    .layer(Bars::new(species, medians))
    .title("Agg::by(...).reduce(Reducer::Median)")
    .y_label("g")

use malevich::stat::ewma;
use malevich::{Line, Plot};
use super::loss_log;
// A real training log, and its debiased exponential moving average.
let loss = loss_log();
Plot::new()
    .layer(Line::y(loss.clone()).label("per-step loss"))
    .layer(Line::y(ewma(&loss, 0.05)).label("ewma α = 0.05").glow())
    .title("topos bigram model, 1,000 steps")
    .x_label("step")
    .y_label("loss")

use malevich::stat::lttb;
use malevich::{Line, Plot, Points};
let x: Vec<f64> = (0..3000).map(|i| f64::from(i) * 0.01).collect();
let y: Vec<f64> = x.iter().map(|x| (x * 2.0).sin() * (x * 0.3).cos() * 3.0 + (x * 17.0).sin() * 0.3).collect();
// LTTB keeps 120 visually representative points. It is explicit and inexact —
// a stat the caller applies — never a default reduction.
let (kept_x, kept_y) = lttb(&x, &y, 120);
Plot::new()
    .layer(Line::xy(x, y).label("3,000 points"))
    .layer(Points::xy(kept_x, kept_y).label("120 kept by LTTB"))
    .title("stat::lttb — an opt-in approximation")

use malevich::stat::Fit;
use super::penguins;
let p = penguins();
let (x, y) = (p.of("Gentoo", &p.flipper), p.of("Gentoo", &p.mass));
let fit = Fit::xy(&x, &y);
// The preset draws the points, the line, and the band; the accumulator answers the numbers.
malevich::trend(x, y)
    .title(format!("Gentoo mass by flipper — slope {:.1} g/mm, R² {:.2}", fit.slope().unwrap_or(0.0), fit.r_squared().unwrap_or(0.0)))
    .x_label("flipper, mm")
    .y_label("g")

use malevich::{ecdf_with, EcdfOptions};
let mut unit = super::noise(13);
let latency: Vec<f64> = (0..600).map(|_| (0.35 * super::gaussian(&mut unit) + 3.0).exp()).collect();
// The empirical CDF, with the Dvoretzky–Kiefer–Wolfowitz 95 % band.
ecdf_with(latency, EcdfOptions::new().band(0.05))
    .expect("valid")
    .title("ECDF with a DKW band")
    .x_label("ms")

use malevich::stat::{auc, roc};
use malevich::{Dash, Line, Plot};
let mut unit = super::noise(17);
// Scores for 400 examples: positives score higher on average, with overlap.
let labels: Vec<bool> = (0..400).map(|i| i % 3 == 0).collect();
let scores: Vec<f64> = labels.iter().map(|&positive| super::gaussian(&mut unit) + if positive { 1.4 } else { 0.0 }).collect();
let (fpr, tpr) = roc(&scores, &labels);
let area = auc(&fpr, &tpr);
Plot::new()
    .layer(Line::xy(fpr, tpr).label("classifier"))
    .layer(Line::xy(vec![0.0, 1.0], vec![0.0, 1.0]).dash(Dash::Dotted).label("chance"))
    .title(format!("ROC, AUC = {area:.3}"))
    .x_label("false positive rate")
    .y_label("true positive rate")

use malevich::stat::stack;
use malevich::{Area, Plot};
let years: Vec<f64> = (2000..2026).map(f64::from).collect();
let coal: Vec<f64> = years.iter().map(|y| 40.0 - (y - 2000.0) * 1.1).collect();
let gas: Vec<f64> = years.iter().map(|y| 20.0 + (y - 2000.0) * 0.4).collect();
let wind: Vec<f64> = years.iter().map(|y| 1.0 + ((y - 2000.0) * 0.35).powi(2)).collect();
let solar: Vec<f64> = years.iter().map(|y| 0.2 + ((y - 2000.0) * 0.28).powi(2)).collect();
// stack returns (low, high) per series; each Area sits on the sum of the ones below.
let mut plot = Plot::new();
for ((low, high), label) in stack(&[&coal, &gas, &wind, &solar]).into_iter().zip(["coal", "gas", "wind", "solar"]) {
    plot = plot.layer(Area::between(years.clone(), low, high).label(label));
}
plot.title("stat::stack — a stacked area").x_label("year").y_label("TWh")

use malevich::stat::{stack_with, StackOffset, StackOptions};
use malevich::{Area, Plot};
let years: Vec<f64> = (2000..2026).map(f64::from).collect();
let coal: Vec<f64> = years.iter().map(|y| 40.0 - (y - 2000.0) * 1.1).collect();
let gas: Vec<f64> = years.iter().map(|y| 20.0 + (y - 2000.0) * 0.4).collect();
let wind: Vec<f64> = years.iter().map(|y| 1.0 + ((y - 2000.0) * 0.35).powi(2)).collect();
let solar: Vec<f64> = years.iter().map(|y| 0.2 + ((y - 2000.0) * 0.28).powi(2)).collect();
// Normalize: every x sums to one, so the chart reads as shares.
let options = StackOptions::new().offset(StackOffset::Normalize);
let mut plot = Plot::new();
for ((low, high), label) in stack_with(&[&coal, &gas, &wind, &solar], options).into_iter().zip(["coal", "gas", "wind", "solar"]) {
    plot = plot.layer(Area::between(years.clone(), low, high).label(label));
}
plot.title("StackOffset::Normalize — shares").x_label("year")

use malevich::stat::dodge;
use malevich::{Bars, Plot, Scale};
let cities = ["Oslo", "Lima", "Osaka", "Cairo"];
let morning = vec![14.0, 22.0, 9.0, 31.0];
let evening = vec![11.0, 19.0, 13.0, 27.0];
// dodge computes side-by-side positions within each band; Bars::at draws there.
let width = 0.38;
let positions = dodge(&[&morning, &evening], width);
Plot::new()
    .x_scale(Scale::bands(cities))
    .layer(Bars::at(positions[0].clone(), width, morning).label("morning"))
    .layer(Bars::at(positions[1].clone(), width, evening).label("evening"))
    .title("stat::dodge — grouped bars")

use malevich::stat::jitter;
use malevich::{Plot, Points, Scale};
use super::penguins;
// A strip plot: every measurement, spread across its band by van der Corput
// offsets — even, seedless, the same every time.
let p = penguins();
let species = ["Adelie", "Chinstrap", "Gentoo"];
let mut plot = Plot::new().x_scale(Scale::bands(species));
for (index, name) in species.iter().enumerate() {
    let masses = p.of(name, &p.mass);
    let positions = jitter(&vec![index as f64; masses.len()], 0.7);
    plot = plot.layer(Points::xy(positions, masses));
}
plot.title("stat::jitter").y_label("g")

use malevich::stat::{steps, StepDirection};
use malevich::{Line, Plot, Points};
let x: Vec<f64> = (0..8).map(f64::from).collect();
let y = vec![1.0, 3.0, 2.0, 5.0, 4.0, 6.0, 3.0, 4.0];
let (ax, ay) = steps(&x, &y, StepDirection::Post);
let (bx, by) = steps(&x, &y, StepDirection::Pre);
let (mx, my) = steps(&x, &y, StepDirection::Mid);
Plot::new()
    .layer(Points::xy(x, y).label("samples"))
    .layer(Line::xy(ax, ay).label("Post"))
    .layer(Line::xy(bx, by.iter().map(|y| y + 6.0).collect::<Vec<_>>()).label("Pre (+6)"))
    .layer(Line::xy(mx, my.iter().map(|y| y + 12.0).collect::<Vec<_>>()).label("Mid (+12)"))
    .title("stat::steps — where a value changes")

use malevich::stat::{cumsum, normalize, Reducer};
use malevich::{Line, Plot};
let mut unit = super::noise(23);
let daily: Vec<f64> = (0..90).map(|i| 3.0 + (f64::from(i) * 0.15).sin() * 2.0 + unit()).collect();
// The series maps: a running total, and a series divided by a reducer of itself.
Plot::new()
    .layer(Line::y(normalize(&daily, Reducer::Max)).label("daily, as a share of the peak"))
    .layer(Line::y(normalize(&cumsum(&daily), Reducer::Max)).label("cumulative, share of the total"))
    .title("cumsum and normalize")
    .x_label("day")

// The MATLAB peaks function as iso-lines: marching squares, tick-chosen levels.
let n = 60;
let mut values = Vec::with_capacity(n * n);
for row in 0..n {
    for column in 0..n {
        let (x, y) = (column as f64 / n as f64 * 6.0 - 3.0, 3.0 - row as f64 / n as f64 * 6.0);
        values.push(3.0 * (1.0 - x).powi(2) * (-x * x - (y + 1.0).powi(2)).exp() - 10.0 * (x / 5.0 - x.powi(3) - y.powi(5)) * (-x * x - y * y).exp() - (-(x + 1.0).powi(2) - y * y).exp() / 3.0);
    }
}
malevich::contour(n, values).title("contour — peaks")

use malevich::{Line, Plot};
let mut unit = super::noise(21);
// 200,000 points with three one-sample spikes. Past four points per column
// the plot reduces by M4 automatically: first, last, min, and max per
// rendered column — the spikes cannot vanish, by construction.
let n = 200_000;
let signal: Vec<f64> = (0..n).map(|i| {
    let t = i as f64 / n as f64;
    let wave = (t * 40.0).sin() * (t * 3.0).cos() * 2.0 + super::gaussian(&mut unit) * 0.2;
    if i == 31_337 || i == 99_999 || i == 150_001 { wave + 7.0 } else { wave }
}).collect();
Plot::new().layer(Line::y(signal)).title("200,000 points through M4").y_min(-4.0).y_max(8.0)

use malevich::{Line, Plot};
let mut unit = super::noise(21);
// The same series, "reduced" the way a sampler would: every 400th point.
// The spikes are gone, and nothing on the chart says so.
let n = 200_000;
let signal: Vec<f64> = (0..n).map(|i| {
    let t = i as f64 / n as f64;
    let wave = (t * 40.0).sin() * (t * 3.0).cos() * 2.0 + super::gaussian(&mut unit) * 0.2;
    if i == 31_337 || i == 99_999 || i == 150_001 { wave + 7.0 } else { wave }
}).collect();
let x: Vec<f64> = (0..n).step_by(400).map(|i| i as f64).collect();
let sampled: Vec<f64> = signal.iter().step_by(400).cloned().collect();
Plot::new().layer(Line::xy(x, sampled)).title("the same points, every 400th").y_min(-4.0).y_max(8.0)