Random distributions

use random_dist.qweruse dist_normal, dist_lognormal from random_dist      // Quillo 0.6: only these

rand = "0.8"rand_distr = "0.4"

// A bakery's Monte Carlo forecast: customers a day (Poisson), what each spends (log-normal),// and the chance a day's takings beat €900; seeded, so the page shows the same run each time.////   quillo run random_dist.qwe      then open http://localhost:8080/forecastuse "random_dist.qwer"route GET /forecast ? { days: int = 2000, seed: int = 42 } {    let customers = dist_poisson(120.0, days, seed)?    let spend = dist_lognormal(1.9, 0.4, days, seed + 1)?    let takings = zip(customers, spend) | map((d) => d[0] * d[1])    let s = dist_summary(takings)?    { typical_day: round(s.p50, 0), bad_day: round(s.p5, 0), good_day: round(s.p95, 0), chance_over_900: round(count(takings, (t) => t > 900) / days, 3) }}route GET /check = {    normal: dist_summary(dist_normal(100.0, 15.0, 100000, 7)?)?,    same_again: dist_normal(0.0, 1.0, 3, 5)? == dist_normal(0.0, 1.0, 3, 5)?,    arrivals_mean: dist_summary(dist_exponential(0.5, 10000, 3)?)?.mean,    defects: dist_binomial(200, 0.02, 5, 9)?,    bad: dist_normal(0.0, -1.0, 3, 1).error,}

//! Random samples from the distributions simulations need: normal, log-normal, Poisson, exponential and binomial, seeded so a run can be repeated, with a summary of what came out. For load tests, test data, Monte Carlo forecasts and games.//! version: 1.0.0//!//! [dependencies]//! rand = "0.8"//! rand_distr = "0.4"use rand::SeedableRng;use rand_distr::Distribution;use serde_json::{json, Value};const MAX_N: i64 = 1_000_000;/// A generator from `seed` (the same seed gives the same numbers on every machine), or from the system's randomness when `seed` is negative.fn rng(seed: i64) -> rand::rngs::StdRng {    if seed < 0 { rand::rngs::StdRng::from_entropy() } else { rand::rngs::StdRng::seed_from_u64(seed as u64) }}fn count(n: i64) -> Result<usize, String> {    if !(1..=MAX_N).contains(&n) {        return Err(format!("n is 1 to {}", MAX_N));    }    Ok(n as usize)}fn take<D: Distribution<f64>>(d: D, n: i64, seed: i64) -> Result<Vec<f64>, String> {    let n = count(n)?;    let mut r = rng(seed);    Ok((0..n).map(|_| d.sample(&mut r)).collect())}/// `n` numbers from a normal (bell curve) distribution with this `mean` and standard deviation `sd` (0 or more). `seed` 0 or more repeats a run; -1 is fresh each time.pub fn dist_normal(mean: f64, sd: f64, n: i64, seed: i64) -> Result<Vec<f64>, String> {    if !(sd >= 0.0 && sd.is_finite() && mean.is_finite()) {        return Err("the mean is a number and sd is 0 or more".into());    }    take(rand_distr::Normal::new(mean, sd).map_err(|e| e.to_string())?, n, seed)}/// `n` numbers from a log-normal distribution: always positive, with a long right tail (incomes, file sizes, response times). `mu` and `sigma` are the mean and sd of the logarithm.pub fn dist_lognormal(mu: f64, sigma: f64, n: i64, seed: i64) -> Result<Vec<f64>, String> {    if !(sigma >= 0.0 && sigma.is_finite() && mu.is_finite()) {        return Err("mu is a number and sigma is 0 or more".into());    }    take(rand_distr::LogNormal::new(mu, sigma).map_err(|e| e.to_string())?, n, seed)}/// `n` counts from a Poisson distribution with mean `lambda` (above 0): how many orders, calls or arrivals in a period that has `lambda` on average.pub fn dist_poisson(lambda: f64, n: i64, seed: i64) -> Result<Vec<i64>, String> {    Ok(take(rand_distr::Poisson::new(lambda).map_err(|e| e.to_string())?, n, seed)?.into_iter().map(|x| x as i64).collect())}/// `n` waiting times from an exponential distribution with `rate` events per unit of time (above 0): the time between arrivals; their mean is 1 / rate.pub fn dist_exponential(rate: f64, n: i64, seed: i64) -> Result<Vec<f64>, String> {    take(rand_distr::Exp::new(rate).map_err(|e| e.to_string())?, n, seed)}/// `n` counts of successes in `trials` tries that each succeed with probability `p` (0 to 1): conversions out of visitors, defects in a batch.pub fn dist_binomial(trials: i64, p: f64, n: i64, seed: i64) -> Result<Vec<i64>, String> {    if trials < 0 {        return Err("trials is 0 or more".into());    }    let d = rand_distr::Binomial::new(trials as u64, p).map_err(|e| e.to_string())?;    let n = count(n)?;    let mut r = rng(seed);    Ok((0..n).map(|_| d.sample(&mut r) as i64).collect())}/// What a list of numbers looks like: `{ n, mean, sd, min, max, p5, p50, p95 }` (percentiles interpolated; sd of the sample).pub fn dist_summary(values: Vec<f64>) -> Result<Value, String> {    let mut v: Vec<f64> = values.into_iter().filter(|x| x.is_finite()).collect();    if v.is_empty() {        return Err("no finite numbers to summarise".into());    }    v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));    let n = v.len() as f64;    let mean = v.iter().sum::<f64>() / n;    let sd = if v.len() > 1 { (v.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / (n - 1.0)).sqrt() } else { 0.0 };    let pct = |p: f64| {        let i = p * (n - 1.0);        let (lo, hi) = (i.floor() as usize, i.ceil() as usize);        v[lo] + (v[hi] - v[lo]) * (i - lo as f64)    };    let r = |x: f64| (x * 1e4).round() / 1e4;    Ok(json!({ "n": v.len(), "mean": r(mean), "sd": r(sd), "min": r(v[0]), "max": r(v[v.len() - 1]), "p5": r(pct(0.05)), "p50": r(pct(0.5)), "p95": r(pct(0.95)) }))}