mirror of
https://github.com/storytold/gridcraft.git
synced 2026-11-03 00:41:56 -06:00
prxnvt
9f4ca54e9a
A Python cell whose cached value is empty (a file's PY formula with no cached value) was taken for an empty cell by six checks, which use the value and the formula only: - the Analysis ToolPak report and the PivotTable report replaced it without asking for overwrite or replace (CodeRabbit, on the report); - Insert Cells and Insert Table without headers pushed it off the end of the sheet; - the data region (Ctrl+A, AutoSum, sort) ended at it; - pasting it into part of a legacy array was not refused. They now use Cell::has_content, which counts Python code. Test: python_code_with_an_empty_value_counts_as_content; reverting any one of the six changes makes it fail.
2590 lines
122 KiB
Rust
2590 lines
122 KiB
Rust
//! Analysis ToolPak (Data › Analysis › Data Analysis): statistical tools that compute a report
|
||
//! from input ranges and write it as a static table (values, not formulas) to an output range or
|
||
//! a new worksheet.
|
||
//!
|
||
//! Every tool is one `tool` value of `data.analysisToolPak`. Distribution values (t, F) come
|
||
//! from the worksheet functions (`T.DIST.2T`, `F.DIST.RT`, `T.INV.2T`…), so the reports agree
|
||
//! with the formulas a user would write. Layouts follow Microsoft's public description of the
|
||
//! tools (row labels and section titles); exact spacing is GridCraft's reading of it.
|
||
|
||
use std::collections::{BTreeMap, HashMap};
|
||
|
||
use gridcraft_core::{CellError, CellRef, DateSystem, MAX_COLS, MAX_ROWS, RangeRef, Value};
|
||
use gridcraft_model::{Cell, NumFmt, Sheet};
|
||
use serde_json::{Value as Json, json};
|
||
|
||
use super::whatif::parse_range;
|
||
use super::{CommandSpec, EngineError, Result, Session, bad, cmd, edit, has_doc, str_param};
|
||
|
||
const CMD: &str = "data.analysisToolPak";
|
||
/// Most input cells a tool reads.
|
||
const MAX_INPUT_CELLS: u64 = 1_000_000;
|
||
/// Most cells a report may write.
|
||
const MAX_OUTPUT_CELLS: u64 = 4_000_000;
|
||
/// Most variables in a correlation or covariance matrix (work grows with its square).
|
||
const MAX_MATRIX_VARS: usize = 256;
|
||
/// Excel's limit on regression X columns.
|
||
const MAX_REGRESSORS: usize = 16;
|
||
/// Most histogram bins.
|
||
const MAX_BINS: usize = 100_000;
|
||
/// Excel's limit on Fourier Analysis input values.
|
||
const MAX_FOURIER: usize = 4096;
|
||
/// Most numbers Random Number Generation or Sampling writes in one run.
|
||
const MAX_RANDOM: usize = 1_000_000;
|
||
/// Most trials of a binomial draw, and the farthest Poisson outcome considered.
|
||
const MAX_OUTCOMES: f64 = 1_000_000.0;
|
||
/// Largest seed: seeds are whole numbers that a JSON number holds exactly.
|
||
const MAX_SEED: u64 = (1 << 53) - 1;
|
||
|
||
/// Part of the error a programmatic call gets when the output range already holds data; the
|
||
/// dialog recognises it and asks before running again with `overwrite: true`.
|
||
pub const OVERWRITE_ERROR: &str = "already holds data";
|
||
|
||
/// (id, label) of every tool, in the order of Excel's Analysis Tools list.
|
||
pub const TOOLS: &[(&str, &str)] = &[
|
||
("anovaSingle", "Anova: Single Factor"),
|
||
("anovaTwoFactorWith", "Anova: Two-Factor With Replication"),
|
||
("anovaTwoFactorWithout", "Anova: Two-Factor Without Replication"),
|
||
("correlation", "Correlation"),
|
||
("covariance", "Covariance"),
|
||
("descriptive", "Descriptive Statistics"),
|
||
("exponentialSmoothing", "Exponential Smoothing"),
|
||
("fTest", "F-Test Two-Sample for Variances"),
|
||
("fourierAnalysis", "Fourier Analysis"),
|
||
("histogram", "Histogram"),
|
||
("movingAverage", "Moving Average"),
|
||
("randomNumberGeneration", "Random Number Generation"),
|
||
("rankPercentile", "Rank and Percentile"),
|
||
("regression", "Regression"),
|
||
("sampling", "Sampling"),
|
||
("tTestPaired", "t-Test: Paired Two Sample for Means"),
|
||
("tTestEqual", "t-Test: Two-Sample Assuming Equal Variances"),
|
||
("tTestUnequal", "t-Test: Two-Sample Assuming Unequal Variances"),
|
||
("zTest", "z-Test: Two Sample for Means"),
|
||
];
|
||
|
||
pub fn specs() -> Vec<CommandSpec> {
|
||
vec![cmd!(
|
||
"data.analysisToolPak",
|
||
"Data Analysis…",
|
||
["Data", "Analysis"],
|
||
None,
|
||
"{tool: anovaSingle|anovaTwoFactorWith|anovaTwoFactorWithout|correlation|covariance|descriptive|exponentialSmoothing|fTest|fourierAnalysis|histogram|movingAverage|randomNumberGeneration|rankPercentile|regression|sampling|tTestPaired|tTestEqual|tTestUnequal|zTest (or the tool's label), output?: \"E1\" (top-left of the report; default a new sheet), newSheet?: name, outputTo?: range|newSheet (picks between the two), overwrite?: bool (else non-empty output cells are refused), labels?: bool (first row/column holds names), groupedBy?: columns|rows, …tool params} — range tools take `range`; descriptive: summary? (default true), confidence?: 95, kthLargest?, kthSmallest?; histogram: bins?: range (default evenly spaced), pareto?, cumulative?; regression: inputY, inputX, constantZero?, confidence?: 95, residuals?, standardizedResiduals?, normalProbability?; movingAverage: interval? (3), standardErrors?; exponentialSmoothing: damping? (0.3), standardErrors?; t-tests and fTest: range1, range2, alpha? (0.05), hypothesizedDifference? (0); anovaSingle: alpha?; anovaTwoFactorWithout: range, labels? (the first row holds column names and the first column row names), alpha?; anovaTwoFactorWith: range (first row and column hold labels), rowsPerSample, alpha?; zTest: range1, range2, variance1, variance2 (known variances), alpha?, hypothesizedDifference?; fourierAnalysis: range (one row or column of 2^k numbers or complex text, at most 4096), inverse?; randomNumberGeneration: variables? (1), count (per variable; not for patterned), distribution: uniform{from, to} | normal{mean, standardDeviation} | bernoulli{p} | binomial{p, trials} | poisson{lambda} | patterned{from, to, step, repeatEach, repeatSequence} | discrete{valueRange: two columns of value and probability}, seed? (whole number, not for patterned; the same seed repeats the numbers; the result reports the seed used), at most 1,000,000 numbers; sampling: range, method?: periodic|random, period (periodic), samples (random, with replacement), seed?. Writes static values. No `tool` → {tools} (the menu opens the dialog).",
|
||
has_doc,
|
||
run
|
||
)]
|
||
}
|
||
|
||
// ---------------------------------------------------------------- parameters
|
||
|
||
fn num(p: &Json, key: &str) -> Option<f64> {
|
||
match p.get(key)? {
|
||
Json::Number(n) => n.as_f64(),
|
||
Json::String(s) => s.trim().trim_end_matches('%').trim().parse::<f64>().ok(),
|
||
_ => None,
|
||
}
|
||
.filter(|v| v.is_finite())
|
||
}
|
||
|
||
fn flag(p: &Json, key: &str) -> bool {
|
||
match p.get(key) {
|
||
Some(Json::Bool(b)) => *b,
|
||
Some(Json::String(s)) => s.eq_ignore_ascii_case("true") || s == "1",
|
||
Some(Json::Number(n)) => n.as_f64().is_some_and(|v| v != 0.0),
|
||
_ => false,
|
||
}
|
||
}
|
||
|
||
fn text<'a>(p: &'a Json, keys: &[&str]) -> Option<&'a str> {
|
||
keys.iter().find_map(|k| str_param(p, k)).map(str::trim).filter(|t| !t.is_empty())
|
||
}
|
||
|
||
fn norm(s: &str) -> String {
|
||
s.chars().filter(char::is_ascii_alphanumeric).collect::<String>().to_ascii_lowercase()
|
||
}
|
||
|
||
/// A tool id from its id or its label (case and punctuation ignored).
|
||
fn tool_id(t: &str) -> Option<&'static str> {
|
||
let n = norm(t);
|
||
TOOLS.iter().find(|(id, label)| norm(id) == n || norm(label) == n).map(|(id, _)| *id)
|
||
}
|
||
|
||
/// Confidence level as a percentage in (0, 100): `95`, `"95%"` or `0.95`.
|
||
fn confidence(p: &Json, key: &str) -> Result<Option<f64>> {
|
||
match num(p, key) {
|
||
None => Ok(None),
|
||
Some(c) => {
|
||
let c = if c > 0.0 && c < 1.0 { c * 100.0 } else { c };
|
||
if c > 0.0 && c < 100.0 { Ok(Some(c)) } else { Err(bad(CMD, "the confidence level must be between 0 and 100")) }
|
||
}
|
||
}
|
||
}
|
||
|
||
fn alpha(p: &Json) -> Result<f64> {
|
||
match opt_num(p, "alpha")? {
|
||
None => Ok(0.05),
|
||
Some(a) if a > 0.0 && a < 1.0 => Ok(a),
|
||
Some(_) => Err(bad(CMD, "alpha must be between 0 and 1")),
|
||
}
|
||
}
|
||
|
||
// ---------------------------------------------------------------- input
|
||
|
||
struct Var {
|
||
name: String,
|
||
vals: Vec<Value>,
|
||
}
|
||
|
||
fn non_numeric() -> EngineError {
|
||
EngineError::Other("Input range contains non-numeric data.".into())
|
||
}
|
||
|
||
/// The cells of a range parameter, after checking its size.
|
||
fn input_range(s: &Session, p: &Json, keys: &[&str]) -> Result<(usize, RangeRef)> {
|
||
let t = text(p, keys).ok_or_else(|| bad(CMD, format!("missing `{}` (the input range)", keys.first().copied().unwrap_or("range"))))?;
|
||
let (si, r) = parse_range(&s.doc()?.wb, t, CMD)?;
|
||
if r.count() > MAX_INPUT_CELLS {
|
||
return Err(bad(CMD, format!("the input range is too large ({MAX_INPUT_CELLS} cells at most)")));
|
||
}
|
||
Ok((si, r))
|
||
}
|
||
|
||
/// Splits a range into variables (columns, or rows when `by_rows`); with `labels` the first
|
||
/// cell of each names it, else `default(i)` does.
|
||
fn variables(s: &Session, si: usize, r: RangeRef, by_rows: bool, labels: bool, default: &dyn Fn(usize) -> String) -> Result<Vec<Var>> {
|
||
let d = s.doc()?;
|
||
let sh = d.wb.sheet(si).ok_or(EngineError::NoDocument)?;
|
||
let (outer, inner) = if by_rows { (r.height(), r.width()) } else { (r.width(), r.height()) };
|
||
let mut out = Vec::with_capacity(outer as usize);
|
||
for i in 0..outer {
|
||
let cell = |j: u32| if by_rows { CellRef::new(r.start.row + i, r.start.col + j) } else { CellRef::new(r.start.row + j, r.start.col + i) };
|
||
let mut vals: Vec<Value> = (0..inner).map(|j| sh.value(cell(j))).collect();
|
||
let mut name = default(i as usize);
|
||
if labels && !vals.is_empty() {
|
||
let first = vals.remove(0);
|
||
let t = first.display();
|
||
if !t.trim().is_empty() {
|
||
name = t;
|
||
}
|
||
}
|
||
out.push(Var { name, vals });
|
||
}
|
||
Ok(out)
|
||
}
|
||
|
||
/// The numbers of a variable; blanks are skipped when `skip_blank`, anything else that isn't a
|
||
/// number is an error.
|
||
fn numbers(v: &[Value], skip_blank: bool) -> Result<Vec<f64>> {
|
||
let mut out = Vec::with_capacity(v.len());
|
||
for x in v {
|
||
match x {
|
||
Value::Number(n) => out.push(*n),
|
||
Value::Empty if skip_blank => {}
|
||
_ => return Err(non_numeric()),
|
||
}
|
||
}
|
||
Ok(out)
|
||
}
|
||
|
||
fn by_rows(p: &Json) -> bool {
|
||
str_param(p, "groupedBy").is_some_and(|g| g.to_ascii_lowercase().starts_with("row"))
|
||
}
|
||
|
||
/// One variable from a single row or column (`range1`/`range2` of the two-sample tests).
|
||
fn single(s: &Session, p: &Json, key: &str, fallback: &str) -> Result<Var> {
|
||
let (si, r) = input_range(s, p, &[key])?;
|
||
if r.width() > 1 && r.height() > 1 {
|
||
return Err(EngineError::Other("Input range must be a single row or column.".into()));
|
||
}
|
||
let rows = r.width() > 1;
|
||
let mut v = variables(s, si, r, rows, flag(p, "labels"), &|_| fallback.to_string())?;
|
||
v.pop().ok_or_else(|| bad(CMD, "empty input range"))
|
||
}
|
||
|
||
// ---------------------------------------------------------------- statistics
|
||
|
||
fn mean(x: &[f64]) -> f64 {
|
||
if x.is_empty() { f64::NAN } else { x.iter().sum::<f64>() / x.len() as f64 }
|
||
}
|
||
|
||
/// Sum of squared deviations from the mean.
|
||
fn ss(x: &[f64]) -> f64 {
|
||
let m = mean(x);
|
||
x.iter().map(|v| (v - m) * (v - m)).sum()
|
||
}
|
||
|
||
/// Sample variance (n − 1).
|
||
fn var_s(x: &[f64]) -> f64 {
|
||
if x.len() < 2 { f64::NAN } else { ss(x) / (x.len() - 1) as f64 }
|
||
}
|
||
|
||
fn div0() -> Value {
|
||
Value::Error(CellError::Div0)
|
||
}
|
||
|
||
/// `a / b`, `#DIV/0!` when `b` is zero or the inputs aren't finite.
|
||
fn ratio(a: f64, b: f64) -> Value {
|
||
if b == 0.0 || !a.is_finite() || !b.is_finite() { div0() } else { Value::number(a / b) }
|
||
}
|
||
|
||
fn n_or_div0(x: f64) -> Value {
|
||
if x.is_finite() { Value::number(x) } else { div0() }
|
||
}
|
||
|
||
/// Host services for calling worksheet functions directly.
|
||
struct FnHost;
|
||
|
||
impl gridcraft_functions::Ctx for FnHost {
|
||
fn date_system(&self) -> DateSystem {
|
||
DateSystem::D1900
|
||
}
|
||
fn now_serial(&self) -> f64 {
|
||
0.0
|
||
}
|
||
fn random(&mut self) -> f64 {
|
||
0.5
|
||
}
|
||
}
|
||
|
||
/// Calls a worksheet function on numeric arguments (and `cumulative` flags as booleans).
|
||
fn call(name: &str, args: &[Value]) -> Value {
|
||
let Some(spec) = gridcraft_functions::lookup(name) else { return Value::Error(CellError::Name) };
|
||
if args.iter().any(|a| a.as_f64().is_some_and(|v| !v.is_finite())) {
|
||
return Value::Error(CellError::Num);
|
||
}
|
||
let a: Vec<gridcraft_functions::Arg> = args.iter().cloned().map(gridcraft_functions::Arg::val).collect();
|
||
gridcraft_functions::call(spec, &a, &mut FnHost)
|
||
}
|
||
|
||
fn callf(name: &str, args: &[f64]) -> Value {
|
||
call(name, &args.iter().map(|v| Value::Number(*v)).collect::<Vec<_>>())
|
||
}
|
||
|
||
fn f64_of(v: &Value) -> f64 {
|
||
v.as_f64().unwrap_or(f64::NAN)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- report grid
|
||
|
||
/// A report: values by (row, col) relative to its top-left corner.
|
||
#[derive(Default)]
|
||
struct Grid {
|
||
cells: BTreeMap<(u32, u32), Value>,
|
||
/// Cells shown as percentages.
|
||
pct: Vec<(u32, u32)>,
|
||
}
|
||
|
||
impl Grid {
|
||
fn put(&mut self, r: u32, c: u32, v: Value) {
|
||
if !v.is_empty() {
|
||
self.cells.insert((r, c), v);
|
||
}
|
||
}
|
||
fn txt(&mut self, r: u32, c: u32, t: impl AsRef<str>) {
|
||
self.put(r, c, Value::text(t.as_ref()));
|
||
}
|
||
fn num(&mut self, r: u32, c: u32, x: f64) {
|
||
self.put(r, c, Value::number(x));
|
||
}
|
||
/// Writes a row of labels starting at column `c`.
|
||
fn labels(&mut self, r: u32, c: u32, items: &[&str]) {
|
||
for (k, t) in items.iter().enumerate() {
|
||
if !t.is_empty() {
|
||
self.txt(r, c + k as u32, t);
|
||
}
|
||
}
|
||
}
|
||
/// (rows, cols) spanned.
|
||
fn extent(&self) -> (u32, u32) {
|
||
let rows = self.cells.keys().map(|(r, _)| r + 1).max().unwrap_or(0);
|
||
let cols = self.cells.keys().map(|(_, c)| c + 1).max().unwrap_or(0);
|
||
(rows, cols)
|
||
}
|
||
}
|
||
|
||
fn idx(i: usize) -> u32 {
|
||
u32::try_from(i).unwrap_or(u32::MAX)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- tools
|
||
|
||
fn descriptive(s: &Session, p: &Json) -> Result<Grid> {
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
let vars = variables(s, si, r, by_rows(p), flag(p, "labels"), &|i| format!("Column{}", i + 1))?;
|
||
let summary = p.get("summary").map(|_| flag(p, "summary")).unwrap_or(true);
|
||
let conf = confidence(p, "confidence")?;
|
||
let kl = num(p, "kthLargest").filter(|k| *k >= 1.0).map(|k| k as usize);
|
||
let ks = num(p, "kthSmallest").filter(|k| *k >= 1.0).map(|k| k as usize);
|
||
if !summary && conf.is_none() && kl.is_none() && ks.is_none() {
|
||
return Err(bad(CMD, "choose at least one of summary, confidence, kthLargest or kthSmallest"));
|
||
}
|
||
let mut g = Grid::default();
|
||
for (i, v) in vars.iter().enumerate() {
|
||
let x = numbers(&v.vals, true)?;
|
||
let c = idx(i).saturating_mul(2);
|
||
g.txt(0, c, &v.name);
|
||
let mut row = 2u32;
|
||
let mut line = |g: &mut Grid, label: &str, val: Value| {
|
||
g.txt(row, c, label);
|
||
g.put(row, c + 1, val);
|
||
row += 1;
|
||
};
|
||
let n = x.len() as f64;
|
||
let mut sorted = x.clone();
|
||
sorted.sort_by(f64::total_cmp);
|
||
if summary {
|
||
let m = mean(&x);
|
||
let sd = var_s(&x).sqrt();
|
||
line(&mut g, "Mean", n_or_div0(m));
|
||
line(&mut g, "Standard Error", if x.len() < 2 { div0() } else { n_or_div0(sd / n.sqrt()) });
|
||
line(&mut g, "Median", median(&sorted));
|
||
line(&mut g, "Mode", mode(&x));
|
||
line(&mut g, "Standard Deviation", n_or_div0(sd));
|
||
line(&mut g, "Sample Variance", n_or_div0(var_s(&x)));
|
||
line(&mut g, "Kurtosis", kurtosis(&x));
|
||
line(&mut g, "Skewness", skewness(&x));
|
||
let (lo, hi) = (sorted.first().copied(), sorted.last().copied());
|
||
line(&mut g, "Range", lo.zip(hi).map(|(a, b)| Value::number(b - a)).unwrap_or(Value::Number(0.0)));
|
||
line(&mut g, "Minimum", lo.map(Value::number).unwrap_or(Value::Number(0.0)));
|
||
line(&mut g, "Maximum", hi.map(Value::number).unwrap_or(Value::Number(0.0)));
|
||
line(&mut g, "Sum", Value::number(x.iter().sum()));
|
||
line(&mut g, "Count", Value::number(n));
|
||
}
|
||
if let Some(k) = kl {
|
||
let v = if k <= sorted.len() { sorted.get(sorted.len() - k).copied().map(Value::number) } else { None };
|
||
line(&mut g, &format!("Largest({k})"), v.unwrap_or(Value::Error(CellError::Num)));
|
||
}
|
||
if let Some(k) = ks {
|
||
let v = sorted.get(k - 1).copied().map(Value::number);
|
||
line(&mut g, &format!("Smallest({k})"), v.unwrap_or(Value::Error(CellError::Num)));
|
||
}
|
||
if let Some(cl) = conf {
|
||
let sd = var_s(&x).sqrt();
|
||
let val = if x.len() < 2 {
|
||
div0()
|
||
} else {
|
||
let t = f64_of(&callf("T.INV.2T", &[1.0 - cl / 100.0, n - 1.0]));
|
||
n_or_div0(t * sd / n.sqrt())
|
||
};
|
||
line(&mut g, &format!("Confidence Level({cl:.1}%)"), val);
|
||
}
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
fn median(sorted: &[f64]) -> Value {
|
||
let n = sorted.len();
|
||
if n == 0 {
|
||
return Value::Error(CellError::Num);
|
||
}
|
||
let v = if n % 2 == 1 { sorted.get(n / 2).copied() } else { sorted.get(n / 2 - 1).zip(sorted.get(n / 2)).map(|(a, b)| (a + b) / 2.0) };
|
||
v.map(Value::number).unwrap_or(Value::Error(CellError::Num))
|
||
}
|
||
|
||
/// The most frequent value (the first one to reach the top count), `#N/A` when no value repeats.
|
||
fn mode(x: &[f64]) -> Value {
|
||
let mut counts: HashMap<u64, usize> = HashMap::new();
|
||
for v in x {
|
||
*counts.entry(v.to_bits()).or_default() += 1;
|
||
}
|
||
let top = counts.values().copied().max().unwrap_or(0);
|
||
if top < 2 {
|
||
return Value::Error(CellError::NA);
|
||
}
|
||
x.iter().find(|v| counts.get(&v.to_bits()).copied() == Some(top)).map(|v| Value::number(*v)).unwrap_or(Value::Error(CellError::NA))
|
||
}
|
||
|
||
/// Sample excess kurtosis (Excel's KURT).
|
||
fn kurtosis(x: &[f64]) -> Value {
|
||
let n = x.len() as f64;
|
||
let sd = var_s(x).sqrt();
|
||
if x.len() < 4 || sd.is_nan() || sd <= 0.0 {
|
||
return div0();
|
||
}
|
||
let m = mean(x);
|
||
let s4: f64 = x.iter().map(|v| ((v - m) / sd).powi(4)).sum();
|
||
Value::number(n * (n + 1.0) / ((n - 1.0) * (n - 2.0) * (n - 3.0)) * s4 - 3.0 * (n - 1.0) * (n - 1.0) / ((n - 2.0) * (n - 3.0)))
|
||
}
|
||
|
||
/// Sample skewness (Excel's SKEW).
|
||
fn skewness(x: &[f64]) -> Value {
|
||
let n = x.len() as f64;
|
||
let sd = var_s(x).sqrt();
|
||
if x.len() < 3 || sd.is_nan() || sd <= 0.0 {
|
||
return div0();
|
||
}
|
||
let m = mean(x);
|
||
let s3: f64 = x.iter().map(|v| ((v - m) / sd).powi(3)).sum();
|
||
Value::number(n / ((n - 1.0) * (n - 2.0)) * s3)
|
||
}
|
||
|
||
fn histogram(s: &Session, p: &Json) -> Result<Grid> {
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
let labels = flag(p, "labels");
|
||
let d = s.doc()?;
|
||
let sh = d.wb.sheet(si).ok_or(EngineError::NoDocument)?;
|
||
let mut cells: Vec<Value> = r.iter().map(|c| sh.value(c)).collect();
|
||
if labels && !cells.is_empty() {
|
||
cells.remove(0);
|
||
}
|
||
let data = numbers(&cells, true)?;
|
||
let mut bins: Vec<f64> = match text(p, &["bins", "binRange"]) {
|
||
Some(_) => {
|
||
let (bsi, br) = input_range(s, p, &["bins", "binRange"])?;
|
||
let bsh = d.wb.sheet(bsi).ok_or(EngineError::NoDocument)?;
|
||
let mut b: Vec<Value> = br.iter().map(|c| bsh.value(c)).collect();
|
||
if labels && !b.is_empty() {
|
||
b.remove(0);
|
||
}
|
||
numbers(&b, true)?
|
||
}
|
||
None => {
|
||
// Evenly spaced bins from the minimum to the maximum, √n intervals.
|
||
let lo = data.iter().copied().fold(f64::INFINITY, f64::min);
|
||
let hi = data.iter().copied().fold(f64::NEG_INFINITY, f64::max);
|
||
if !lo.is_finite() || !hi.is_finite() {
|
||
return Err(EngineError::Other("The input range holds no numbers.".into()));
|
||
}
|
||
let k = ((data.len() as f64).sqrt().floor() as usize).clamp(1, 1000);
|
||
if hi == lo { vec![lo] } else { (0..=k).map(|i| lo + (hi - lo) * i as f64 / k as f64).collect() }
|
||
}
|
||
};
|
||
if bins.is_empty() {
|
||
return Err(EngineError::Other("The bin range holds no numbers.".into()));
|
||
}
|
||
if bins.len() > MAX_BINS {
|
||
return Err(bad(CMD, format!("too many bins ({MAX_BINS} at most)")));
|
||
}
|
||
bins.sort_by(f64::total_cmp);
|
||
let mut freq = vec![0u64; bins.len() + 1];
|
||
for x in &data {
|
||
let i = bins.partition_point(|b| b < x);
|
||
if let Some(f) = freq.get_mut(i) {
|
||
*f += 1;
|
||
}
|
||
}
|
||
let rows: Vec<(Value, u64)> = bins.iter().map(|b| Value::number(*b)).chain(std::iter::once(Value::text("More"))).zip(freq).collect();
|
||
let total = data.len() as f64;
|
||
let cumulative = flag(p, "cumulative");
|
||
let mut g = Grid::default();
|
||
let write = |g: &mut Grid, c0: u32, rows: &[(Value, u64)]| {
|
||
g.labels(0, c0, &["Bin", "Frequency"]);
|
||
if cumulative {
|
||
g.txt(0, c0 + 2, "Cumulative %");
|
||
}
|
||
let mut acc = 0u64;
|
||
for (i, (b, f)) in rows.iter().enumerate() {
|
||
let r = idx(i) + 1;
|
||
g.put(r, c0, b.clone());
|
||
g.num(r, c0 + 1, *f as f64);
|
||
if cumulative {
|
||
acc += f;
|
||
g.put(r, c0 + 2, ratio(acc as f64, total));
|
||
g.pct.push((r, c0 + 2));
|
||
}
|
||
}
|
||
};
|
||
write(&mut g, 0, &rows);
|
||
if flag(p, "pareto") {
|
||
let mut sorted = rows.clone();
|
||
sorted.sort_by_key(|r| std::cmp::Reverse(r.1));
|
||
write(&mut g, if cumulative { 3 } else { 2 }, &sorted);
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
/// Variables of a range whose columns (rows) must all be numbers, for the matrix tools.
|
||
fn matrix_vars(s: &Session, p: &Json) -> Result<Vec<(String, Vec<f64>)>> {
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
let vars = variables(s, si, r, by_rows(p), flag(p, "labels"), &|i| format!("Column {}", i + 1))?;
|
||
if vars.len() > MAX_MATRIX_VARS {
|
||
return Err(bad(CMD, format!("too many variables ({MAX_MATRIX_VARS} at most)")));
|
||
}
|
||
vars.into_iter().map(|v| Ok((v.name, numbers(&v.vals, false)?))).collect()
|
||
}
|
||
|
||
fn matrix(s: &Session, p: &Json, f: impl Fn(&[f64], &[f64]) -> Value, diag: impl Fn(&[f64]) -> Value) -> Result<Grid> {
|
||
let vars = matrix_vars(s, p)?;
|
||
let mut g = Grid::default();
|
||
for (i, (name, x)) in vars.iter().enumerate() {
|
||
let k = idx(i) + 1;
|
||
g.txt(0, k, name);
|
||
g.txt(k, 0, name);
|
||
for (j, (_, y)) in vars.iter().enumerate().take(i + 1) {
|
||
g.put(k, idx(j) + 1, if i == j { diag(x) } else { f(x, y) });
|
||
}
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
/// Sum of cross products of deviations.
|
||
fn sxy(x: &[f64], y: &[f64]) -> f64 {
|
||
let (mx, my) = (mean(x), mean(y));
|
||
x.iter().zip(y).map(|(a, b)| (a - mx) * (b - my)).sum()
|
||
}
|
||
|
||
fn correlation(x: &[f64], y: &[f64]) -> Value {
|
||
let d = (ss(x) * ss(y)).sqrt();
|
||
if x.is_empty() || d.is_nan() || d <= 0.0 { div0() } else { Value::number(sxy(x, y) / d) }
|
||
}
|
||
|
||
fn covariance_p(x: &[f64], y: &[f64]) -> Value {
|
||
if x.is_empty() { div0() } else { Value::number(sxy(x, y) / x.len() as f64) }
|
||
}
|
||
|
||
// ---------------------------------------------------------------- regression
|
||
|
||
/// Inverts a small symmetric matrix by Gauss–Jordan elimination with partial pivoting; `None`
|
||
/// when it is singular (collinear inputs).
|
||
fn invert(mut a: Vec<Vec<f64>>) -> Option<Vec<Vec<f64>>> {
|
||
let n = a.len();
|
||
let scale = (0..n).filter_map(|i| a.get(i).and_then(|r| r.get(i))).fold(0.0f64, |m, v| m.max(v.abs()));
|
||
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 col in 0..n {
|
||
let piv = (col..n).max_by(|&x, &y| {
|
||
let ax = a.get(x).and_then(|r| r.get(col)).map_or(0.0, |v| v.abs());
|
||
let ay = a.get(y).and_then(|r| r.get(col)).map_or(0.0, |v| v.abs());
|
||
ax.total_cmp(&ay)
|
||
})?;
|
||
let pv = *a.get(piv)?.get(col)?;
|
||
if pv.is_nan() || pv.abs() <= 1e-12 * scale.max(f64::MIN_POSITIVE) {
|
||
return None;
|
||
}
|
||
a.swap(col, piv);
|
||
inv.swap(col, piv);
|
||
for j in 0..n {
|
||
*a.get_mut(col)?.get_mut(j)? /= pv;
|
||
*inv.get_mut(col)?.get_mut(j)? /= pv;
|
||
}
|
||
let prow = a.get(col)?.clone();
|
||
let pinv = inv.get(col)?.clone();
|
||
for i in 0..n {
|
||
if i == col {
|
||
continue;
|
||
}
|
||
let f = *a.get(i)?.get(col)?;
|
||
if f == 0.0 {
|
||
continue;
|
||
}
|
||
for j in 0..n {
|
||
*a.get_mut(i)?.get_mut(j)? -= f * prow.get(j)?;
|
||
*inv.get_mut(i)?.get_mut(j)? -= f * pinv.get(j)?;
|
||
}
|
||
}
|
||
}
|
||
Some(inv)
|
||
}
|
||
|
||
struct Fit {
|
||
/// Intercept first, then one coefficient per X column.
|
||
coef: Vec<f64>,
|
||
/// Standard errors in the same order (NaN for a forced-zero intercept).
|
||
se: Vec<f64>,
|
||
predicted: Vec<f64>,
|
||
sse: f64,
|
||
ssr: f64,
|
||
dfr: f64,
|
||
dfe: f64,
|
||
/// Total degrees of freedom (n − 1, or n through the origin).
|
||
dft: f64,
|
||
}
|
||
|
||
/// Ordinary least squares of `y` on the columns `xs`, with an intercept unless `zero`.
|
||
fn ols(y: &[f64], xs: &[Vec<f64>], zero: bool) -> Option<Fit> {
|
||
let n = y.len();
|
||
let k = xs.len();
|
||
let nf = n as f64;
|
||
let (ym, xm): (f64, Vec<f64>) = if zero { (0.0, vec![0.0; k]) } else { (mean(y), xs.iter().map(|x| mean(x)).collect()) };
|
||
let dev = |j: usize, i: usize| xs.get(j).and_then(|x| x.get(i)).copied().unwrap_or(0.0) - xm.get(j).copied().unwrap_or(0.0);
|
||
let mut sxx = vec![vec![0.0; k]; k];
|
||
let mut sxy = vec![0.0; k];
|
||
for i in 0..n {
|
||
let yi = y.get(i).copied().unwrap_or(0.0) - ym;
|
||
for a in 0..k {
|
||
let da = dev(a, i);
|
||
if let Some(v) = sxy.get_mut(a) {
|
||
*v += da * yi;
|
||
}
|
||
for b in 0..k {
|
||
if let Some(v) = sxx.get_mut(a).and_then(|r| r.get_mut(b)) {
|
||
*v += da * dev(b, i);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let inv = invert(sxx)?;
|
||
let b: Vec<f64> = (0..k)
|
||
.map(|a| (0..k).map(|c| inv.get(a).and_then(|r| r.get(c)).copied().unwrap_or(0.0) * sxy.get(c).copied().unwrap_or(0.0)).sum())
|
||
.collect();
|
||
let b0 = if zero { 0.0 } else { ym - b.iter().zip(&xm).map(|(bj, mj)| bj * mj).sum::<f64>() };
|
||
let predicted: Vec<f64> = (0..n)
|
||
.map(|i| b0 + (0..k).map(|j| b.get(j).copied().unwrap_or(0.0) * xs.get(j).and_then(|x| x.get(i)).copied().unwrap_or(0.0)).sum::<f64>())
|
||
.collect();
|
||
let sse: f64 = y.iter().zip(&predicted).map(|(a, p)| (a - p) * (a - p)).sum();
|
||
let ssr: f64 = predicted.iter().map(|p| (p - ym) * (p - ym)).sum();
|
||
let dfr = k as f64;
|
||
let dfe = if zero { nf - dfr } else { nf - dfr - 1.0 };
|
||
let dft = if zero { nf } else { nf - 1.0 };
|
||
let s2 = sse / dfe;
|
||
let mut se = Vec::with_capacity(k + 1);
|
||
if zero {
|
||
se.push(f64::NAN);
|
||
} else {
|
||
let q: f64 = (0..k)
|
||
.map(|a| {
|
||
(0..k)
|
||
.map(|c| {
|
||
xm.get(a).copied().unwrap_or(0.0)
|
||
* inv.get(a).and_then(|r| r.get(c)).copied().unwrap_or(0.0)
|
||
* xm.get(c).copied().unwrap_or(0.0)
|
||
})
|
||
.sum::<f64>()
|
||
})
|
||
.sum();
|
||
se.push((s2 * (1.0 / nf + q)).sqrt());
|
||
}
|
||
for a in 0..k {
|
||
se.push((s2 * inv.get(a).and_then(|r| r.get(a)).copied().unwrap_or(f64::NAN)).sqrt());
|
||
}
|
||
let mut coef = vec![b0];
|
||
coef.extend(b);
|
||
Some(Fit { coef, se, predicted, sse, ssr, dfr, dfe, dft })
|
||
}
|
||
|
||
fn regression(s: &Session, p: &Json) -> Result<Grid> {
|
||
let labels = flag(p, "labels");
|
||
let (ysi, yr) = input_range(s, p, &["inputY", "y"])?;
|
||
let (xsi, xr) = input_range(s, p, &["inputX", "x"])?;
|
||
if yr.width() != 1 {
|
||
return Err(EngineError::Other("Input Y range must be a single column.".into()));
|
||
}
|
||
if xr.height() != yr.height() {
|
||
return Err(EngineError::Other("Input Y range and input X range must have the same number of rows.".into()));
|
||
}
|
||
if xr.width() as usize > MAX_REGRESSORS {
|
||
return Err(EngineError::Other(format!("Input X range must have {MAX_REGRESSORS} columns or fewer.")));
|
||
}
|
||
let yv = variables(s, ysi, yr, false, labels, &|_| "Y".into())?;
|
||
let xv = variables(s, xsi, xr, false, labels, &|i| format!("X Variable {}", i + 1))?;
|
||
let yvar = yv.first().ok_or_else(|| bad(CMD, "empty input Y range"))?;
|
||
let y = numbers(&yvar.vals, false)?;
|
||
let xs: Vec<Vec<f64>> = xv.iter().map(|v| numbers(&v.vals, false)).collect::<Result<_>>()?;
|
||
let zero = flag(p, "constantZero");
|
||
let k = xs.len();
|
||
let n = y.len();
|
||
let params = if zero { k } else { k + 1 };
|
||
if n <= params {
|
||
return Err(EngineError::Other("Too few observations for this many X variables.".into()));
|
||
}
|
||
let fit = ols(&y, &xs, zero).ok_or_else(|| EngineError::Other("The X variables are collinear; the regression can't be computed.".into()))?;
|
||
let conf = confidence(p, "confidence")?.unwrap_or(95.0);
|
||
let sst = fit.ssr + fit.sse;
|
||
let r2 = if sst > 0.0 { fit.ssr / sst } else { f64::NAN };
|
||
let mut g = Grid::default();
|
||
g.txt(0, 0, "SUMMARY OUTPUT");
|
||
g.txt(2, 0, "Regression Statistics");
|
||
let stats: [(&str, Value); 5] = [
|
||
("Multiple R", n_or_div0(r2.sqrt())),
|
||
("R Square", n_or_div0(r2)),
|
||
("Adjusted R Square", n_or_div0(1.0 - (fit.sse / fit.dfe) / (sst / fit.dft))),
|
||
("Standard Error", n_or_div0((fit.sse / fit.dfe).sqrt())),
|
||
("Observations", Value::number(n as f64)),
|
||
];
|
||
for (i, (label, v)) in stats.into_iter().enumerate() {
|
||
g.txt(3 + idx(i), 0, label);
|
||
g.put(3 + idx(i), 1, v);
|
||
}
|
||
g.txt(9, 0, "ANOVA");
|
||
g.labels(10, 1, &["df", "SS", "MS", "F", "Significance F"]);
|
||
let (msr, mse) = (fit.ssr / fit.dfr, fit.sse / fit.dfe);
|
||
let f = msr / mse;
|
||
g.txt(11, 0, "Regression");
|
||
g.num(11, 1, fit.dfr);
|
||
g.num(11, 2, fit.ssr);
|
||
g.num(11, 3, msr);
|
||
g.put(11, 4, n_or_div0(f));
|
||
g.put(11, 5, if f.is_finite() { callf("F.DIST.RT", &[f, fit.dfr, fit.dfe]) } else { div0() });
|
||
g.txt(12, 0, "Residual");
|
||
g.num(12, 1, fit.dfe);
|
||
g.num(12, 2, fit.sse);
|
||
g.num(12, 3, mse);
|
||
g.txt(13, 0, "Total");
|
||
g.num(13, 1, fit.dft);
|
||
g.num(13, 2, sst);
|
||
let lo2 = format!("Lower {conf:.1}%");
|
||
let hi2 = format!("Upper {conf:.1}%");
|
||
g.labels(15, 1, &["Coefficients", "Standard Error", "t Stat", "P-value", "Lower 95%", "Upper 95%", &lo2, &hi2]);
|
||
let t95 = f64_of(&callf("T.INV.2T", &[0.05, fit.dfe]));
|
||
let tc = f64_of(&callf("T.INV.2T", &[1.0 - conf / 100.0, fit.dfe]));
|
||
for (j, (b, se)) in fit.coef.iter().zip(&fit.se).enumerate() {
|
||
let row = 16 + idx(j);
|
||
g.txt(row, 0, if j == 0 { "Intercept" } else { xv.get(j - 1).map_or("", |v| v.name.as_str()) });
|
||
g.num(row, 1, *b);
|
||
if j == 0 && zero {
|
||
for c in 2..=8 {
|
||
g.put(row, c, Value::Error(CellError::NA));
|
||
}
|
||
continue;
|
||
}
|
||
let t = b / se;
|
||
g.put(row, 2, n_or_div0(*se));
|
||
g.put(row, 3, n_or_div0(t));
|
||
g.put(row, 4, if t.is_finite() { callf("T.DIST.2T", &[t.abs(), fit.dfe]) } else { div0() });
|
||
g.num(row, 5, b - t95 * se);
|
||
g.num(row, 6, b + t95 * se);
|
||
g.num(row, 7, b - tc * se);
|
||
g.num(row, 8, b + tc * se);
|
||
}
|
||
let mut row = 16 + idx(fit.coef.len()) + 3;
|
||
let residuals = flag(p, "residuals");
|
||
let standardized = flag(p, "standardizedResiduals");
|
||
if residuals || standardized {
|
||
g.txt(row, 0, "RESIDUAL OUTPUT");
|
||
let pred = format!("Predicted {}", yvar.name);
|
||
let mut head = vec!["Observation", pred.as_str(), "Residuals"];
|
||
if standardized {
|
||
head.push("Standard Residuals");
|
||
}
|
||
g.labels(row + 2, 0, &head);
|
||
// Residuals scaled by their standard deviation (n − 1).
|
||
let sd = (fit.sse / (n as f64 - 1.0)).sqrt();
|
||
for (i, (yi, pi)) in y.iter().zip(&fit.predicted).enumerate() {
|
||
let r = row + 3 + idx(i);
|
||
g.num(r, 0, (i + 1) as f64);
|
||
g.num(r, 1, *pi);
|
||
g.num(r, 2, yi - pi);
|
||
if standardized {
|
||
g.put(r, 3, ratio(yi - pi, sd));
|
||
}
|
||
}
|
||
}
|
||
if flag(p, "normalProbability") {
|
||
let c = if residuals || standardized { if standardized { 5 } else { 4 } } else { 0 };
|
||
if !(residuals || standardized) {
|
||
row = 16 + idx(fit.coef.len()) + 3;
|
||
}
|
||
g.txt(row, c, "PROBABILITY OUTPUT");
|
||
g.labels(row + 2, c, &["Percentile", &yvar.name]);
|
||
let mut sorted = y.clone();
|
||
sorted.sort_by(f64::total_cmp);
|
||
for (i, v) in sorted.iter().enumerate() {
|
||
let r = row + 3 + idx(i);
|
||
g.num(r, c, 100.0 * (i as f64 + 0.5) / n as f64);
|
||
g.num(r, c + 1, *v);
|
||
}
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- smoothing
|
||
|
||
fn series(s: &Session, p: &Json) -> Result<Vec<f64>> {
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
if r.width() > 1 && r.height() > 1 {
|
||
return Err(EngineError::Other("Input range must be a single row or column.".into()));
|
||
}
|
||
let v = variables(s, si, r, r.width() > 1, flag(p, "labels"), &|_| String::new())?;
|
||
numbers(&v.first().map(|v| v.vals.clone()).unwrap_or_default(), false)
|
||
}
|
||
|
||
/// Writes the forecasts (and with `standardErrors` the standard errors over `se_window`
|
||
/// points) down two columns.
|
||
fn forecast_grid(x: &[f64], f: &[Option<f64>], se_window: usize, p: &Json) -> Grid {
|
||
let mut g = Grid::default();
|
||
let na = || Value::Error(CellError::NA);
|
||
for (i, fi) in f.iter().enumerate() {
|
||
g.put(idx(i), 0, fi.map(Value::number).unwrap_or_else(na));
|
||
}
|
||
if flag(p, "standardErrors") {
|
||
// A running sum of squared residuals over the last `se_window` points: O(n), not
|
||
// O(n·window). A point without a forecast restarts it, so every window touching it is #N/A.
|
||
let sq = |j: usize| Some((x.get(j)? - (*f.get(j)?)?).powi(2));
|
||
let (mut acc, mut run) = (0.0, 0usize);
|
||
for i in 0..f.len() {
|
||
match sq(i) {
|
||
Some(t) => {
|
||
acc += t;
|
||
run += 1;
|
||
}
|
||
None => (acc, run) = (0.0, 0),
|
||
}
|
||
if run > se_window {
|
||
acc -= i.checked_sub(se_window).and_then(sq).unwrap_or(0.0);
|
||
run = se_window;
|
||
}
|
||
// Rounding in the subtraction can leave a tiny negative sum.
|
||
let v = (run == se_window).then(|| (acc.max(0.0) / se_window as f64).sqrt());
|
||
g.put(idx(i), 1, v.map(Value::number).unwrap_or_else(na));
|
||
}
|
||
}
|
||
g
|
||
}
|
||
|
||
fn moving_average(s: &Session, p: &Json) -> Result<Grid> {
|
||
let x = series(s, p)?;
|
||
let interval = num(p, "interval").unwrap_or(3.0);
|
||
if interval < 1.0 || interval.fract() != 0.0 || interval > x.len() as f64 {
|
||
return Err(EngineError::Other("The interval must be a whole number from 1 to the number of data points.".into()));
|
||
}
|
||
let k = interval as usize;
|
||
let mut f = Vec::with_capacity(x.len());
|
||
let mut acc = 0.0;
|
||
for (i, v) in x.iter().enumerate() {
|
||
acc += v;
|
||
if let Some(old) = i.checked_sub(k).and_then(|j| x.get(j)) {
|
||
acc -= old;
|
||
}
|
||
f.push(if i + 1 >= k { Some(acc / k as f64) } else { None });
|
||
}
|
||
Ok(forecast_grid(&x, &f, k, p))
|
||
}
|
||
|
||
fn exponential_smoothing(s: &Session, p: &Json) -> Result<Grid> {
|
||
let x = series(s, p)?;
|
||
let damp = num(p, "damping").unwrap_or(0.3);
|
||
if !(0.0..=1.0).contains(&damp) {
|
||
return Err(EngineError::Other("The damping factor must be between 0 and 1.".into()));
|
||
}
|
||
// F(1) = A(0); F(t) = (1 − damping)·A(t−1) + damping·F(t−1).
|
||
let mut f: Vec<Option<f64>> = Vec::with_capacity(x.len());
|
||
for i in 0..x.len() {
|
||
let v = match i {
|
||
0 => None,
|
||
1 => x.first().copied(),
|
||
_ => x.get(i - 1).zip(f.get(i - 1).copied().flatten()).map(|(a, prev)| (1.0 - damp) * a + damp * prev),
|
||
};
|
||
f.push(v);
|
||
}
|
||
Ok(forecast_grid(&x, &f, 3, p))
|
||
}
|
||
|
||
// ---------------------------------------------------------------- rank and percentile
|
||
|
||
fn rank_percentile(s: &Session, p: &Json) -> Result<Grid> {
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
let vars = variables(s, si, r, by_rows(p), flag(p, "labels"), &|i| format!("Column{}", i + 1))?;
|
||
let mut g = Grid::default();
|
||
for (vi, v) in vars.iter().enumerate() {
|
||
let c = idx(vi).saturating_mul(4);
|
||
g.labels(0, c, &["Point", &v.name, "Rank", "Percent"]);
|
||
let mut pts: Vec<(usize, f64)> = Vec::new();
|
||
for (i, val) in v.vals.iter().enumerate() {
|
||
match val {
|
||
Value::Number(n) => pts.push((i + 1, *n)),
|
||
Value::Empty => {}
|
||
_ => return Err(non_numeric()),
|
||
}
|
||
}
|
||
let mut asc: Vec<f64> = pts.iter().map(|p| p.1).collect();
|
||
asc.sort_by(f64::total_cmp);
|
||
let n = asc.len();
|
||
pts.sort_by(|a, b| b.1.total_cmp(&a.1).then(a.0.cmp(&b.0)));
|
||
for (k, (point, x)) in pts.iter().enumerate() {
|
||
let row = idx(k) + 1;
|
||
let below = asc.partition_point(|v| v < x);
|
||
let above = n - asc.partition_point(|v| v <= x);
|
||
g.num(row, c, *point as f64);
|
||
g.num(row, c + 1, *x);
|
||
g.num(row, c + 2, (above + 1) as f64);
|
||
// PERCENTRANK.INC: share of values below, truncated to three digits.
|
||
let pr = if n <= 1 { 1.0 } else { (below as f64 / (n - 1) as f64 * 1000.0 + 1e-9).floor() / 1000.0 };
|
||
g.num(row, c + 3, pr);
|
||
g.pct.push((row, c + 3));
|
||
}
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- tests of means and variances
|
||
|
||
fn two_samples(s: &Session, p: &Json, skip_blank: bool) -> Result<(Var, Var, Vec<f64>, Vec<f64>)> {
|
||
let a = single(s, p, "range1", "Variable 1")?;
|
||
let b = single(s, p, "range2", "Variable 2")?;
|
||
let x = numbers(&a.vals, skip_blank)?;
|
||
let y = numbers(&b.vals, skip_blank)?;
|
||
if x.len() < 2 || y.len() < 2 {
|
||
return Err(EngineError::Other("Each input range needs at least two numbers.".into()));
|
||
}
|
||
Ok((a, b, x, y))
|
||
}
|
||
|
||
/// Title, blank row and the "Variable 1 / Variable 2" header; returns the first data row.
|
||
fn test_head(g: &mut Grid, title: &str, a: &str, b: &str) -> u32 {
|
||
g.txt(0, 0, title);
|
||
g.txt(2, 1, a);
|
||
g.txt(2, 2, b);
|
||
3
|
||
}
|
||
|
||
/// The t statistic rows shared by the three t-tests.
|
||
fn t_rows(g: &mut Grid, mut row: u32, t: f64, df: f64, alpha: f64) {
|
||
let ta = t.abs();
|
||
let lines: [(&str, Value); 5] = [
|
||
("t Stat", n_or_div0(t)),
|
||
("P(T<=t) one-tail", if t.is_finite() { callf("T.DIST.RT", &[ta, df]) } else { div0() }),
|
||
("t Critical one-tail", callf("T.INV", &[1.0 - alpha, df])),
|
||
("P(T<=t) two-tail", if t.is_finite() { callf("T.DIST.2T", &[ta, df]) } else { div0() }),
|
||
("t Critical two-tail", callf("T.INV.2T", &[alpha, df])),
|
||
];
|
||
for (label, v) in lines {
|
||
g.txt(row, 0, label);
|
||
g.put(row, 1, v);
|
||
row += 1;
|
||
}
|
||
}
|
||
|
||
fn pair_row(g: &mut Grid, row: u32, label: &str, a: Value, b: Value) {
|
||
g.txt(row, 0, label);
|
||
g.put(row, 1, a);
|
||
g.put(row, 2, b);
|
||
}
|
||
|
||
fn t_test(s: &Session, p: &Json, kind: &str) -> Result<Grid> {
|
||
let alpha = alpha(p)?;
|
||
let hmd = num(p, "hypothesizedDifference").or_else(|| num(p, "difference")).unwrap_or(0.0);
|
||
let paired = kind == "tTestPaired";
|
||
let (a, b, x, y) = two_samples(s, p, !paired)?;
|
||
let (n1, n2) = (x.len() as f64, y.len() as f64);
|
||
let (m1, m2, v1, v2) = (mean(&x), mean(&y), var_s(&x), var_s(&y));
|
||
let mut g = Grid::default();
|
||
let title = match kind {
|
||
"tTestPaired" => "t-Test: Paired Two Sample for Means",
|
||
"tTestEqual" => "t-Test: Two-Sample Assuming Equal Variances",
|
||
_ => "t-Test: Two-Sample Assuming Unequal Variances",
|
||
};
|
||
let mut row = test_head(&mut g, title, &a.name, &b.name);
|
||
pair_row(&mut g, row, "Mean", Value::number(m1), Value::number(m2));
|
||
pair_row(&mut g, row + 1, "Variance", Value::number(v1), Value::number(v2));
|
||
pair_row(&mut g, row + 2, "Observations", Value::number(n1), Value::number(n2));
|
||
row += 3;
|
||
let (t, df) = match kind {
|
||
"tTestPaired" => {
|
||
if x.len() != y.len() {
|
||
return Err(EngineError::Other("Input ranges must have the same number of values.".into()));
|
||
}
|
||
pair_row(&mut g, row, "Pearson Correlation", correlation(&x, &y), Value::Empty);
|
||
row += 1;
|
||
let dif: Vec<f64> = x.iter().zip(&y).map(|(a, b)| a - b).collect();
|
||
let sd = var_s(&dif).sqrt();
|
||
((mean(&dif) - hmd) / (sd / n1.sqrt()), n1 - 1.0)
|
||
}
|
||
"tTestEqual" => {
|
||
let pooled = ((n1 - 1.0) * v1 + (n2 - 1.0) * v2) / (n1 + n2 - 2.0);
|
||
pair_row(&mut g, row, "Pooled Variance", Value::number(pooled), Value::Empty);
|
||
row += 1;
|
||
((m1 - m2 - hmd) / (pooled * (1.0 / n1 + 1.0 / n2)).sqrt(), n1 + n2 - 2.0)
|
||
}
|
||
_ => {
|
||
let (q1, q2) = (v1 / n1, v2 / n2);
|
||
// Welch–Satterthwaite degrees of freedom, rounded to a whole number.
|
||
let df = ((q1 + q2).powi(2) / (q1 * q1 / (n1 - 1.0) + q2 * q2 / (n2 - 1.0))).round();
|
||
((m1 - m2 - hmd) / (q1 + q2).sqrt(), df)
|
||
}
|
||
};
|
||
pair_row(&mut g, row, "Hypothesized Mean Difference", Value::number(hmd), Value::Empty);
|
||
pair_row(&mut g, row + 1, "df", n_or_div0(df), Value::Empty);
|
||
t_rows(&mut g, row + 2, t, df, alpha);
|
||
Ok(g)
|
||
}
|
||
|
||
fn f_test(s: &Session, p: &Json) -> Result<Grid> {
|
||
let alpha = alpha(p)?;
|
||
let (a, b, x, y) = two_samples(s, p, true)?;
|
||
let (v1, v2) = (var_s(&x), var_s(&y));
|
||
let (d1, d2) = (x.len() as f64 - 1.0, y.len() as f64 - 1.0);
|
||
let mut g = Grid::default();
|
||
let row = test_head(&mut g, "F-Test Two-Sample for Variances", &a.name, &b.name);
|
||
pair_row(&mut g, row, "Mean", Value::number(mean(&x)), Value::number(mean(&y)));
|
||
pair_row(&mut g, row + 1, "Variance", Value::number(v1), Value::number(v2));
|
||
pair_row(&mut g, row + 2, "Observations", Value::number(d1 + 1.0), Value::number(d2 + 1.0));
|
||
pair_row(&mut g, row + 3, "df", Value::number(d1), Value::number(d2));
|
||
let f = v1 / v2;
|
||
pair_row(&mut g, row + 4, "F", n_or_div0(f), Value::Empty);
|
||
// The tail on the side of F: left when F < 1, right otherwise.
|
||
let (pv, crit) = if !f.is_finite() {
|
||
(div0(), div0())
|
||
} else if f < 1.0 {
|
||
(call("F.DIST", &[Value::Number(f), Value::Number(d1), Value::Number(d2), Value::Bool(true)]), callf("F.INV", &[alpha, d1, d2]))
|
||
} else {
|
||
(callf("F.DIST.RT", &[f, d1, d2]), callf("F.INV.RT", &[alpha, d1, d2]))
|
||
};
|
||
pair_row(&mut g, row + 5, "P(F<=f) one-tail", pv, Value::Empty);
|
||
pair_row(&mut g, row + 6, "F Critical one-tail", crit, Value::Empty);
|
||
Ok(g)
|
||
}
|
||
|
||
fn anova_single(s: &Session, p: &Json) -> Result<Grid> {
|
||
let alpha = alpha(p)?;
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
let vars = variables(s, si, r, by_rows(p), flag(p, "labels"), &|i| format!("Column {}", i + 1))?;
|
||
let groups: Vec<(String, Vec<f64>)> = vars.into_iter().map(|v| Ok((v.name, numbers(&v.vals, true)?))).collect::<Result<_>>()?;
|
||
if groups.len() < 2 || groups.iter().any(|(_, x)| x.is_empty()) {
|
||
return Err(EngineError::Other("Anova needs at least two groups, each with at least one number.".into()));
|
||
}
|
||
let mut g = Grid::default();
|
||
g.txt(0, 0, "Anova: Single Factor");
|
||
g.txt(2, 0, "SUMMARY");
|
||
g.labels(3, 0, &["Groups", "Count", "Sum", "Average", "Variance"]);
|
||
let all: Vec<f64> = groups.iter().flat_map(|(_, x)| x.iter().copied()).collect();
|
||
let grand = mean(&all);
|
||
let (mut ssb, mut ssw) = (0.0, 0.0);
|
||
for (i, (name, x)) in groups.iter().enumerate() {
|
||
let row = 4 + idx(i);
|
||
let m = mean(x);
|
||
g.txt(row, 0, name);
|
||
g.num(row, 1, x.len() as f64);
|
||
g.num(row, 2, x.iter().sum());
|
||
g.num(row, 3, m);
|
||
g.put(row, 4, n_or_div0(var_s(x)));
|
||
ssb += x.len() as f64 * (m - grand) * (m - grand);
|
||
ssw += ss(x);
|
||
}
|
||
let k = groups.len() as f64;
|
||
let n = all.len() as f64;
|
||
let (dfb, dfw) = (k - 1.0, n - k);
|
||
let row = 4 + idx(groups.len()) + 2;
|
||
g.txt(row, 0, "ANOVA");
|
||
g.labels(row + 1, 0, &["Source of Variation", "SS", "df", "MS", "F", "P-value", "F crit"]);
|
||
let (msb, msw) = (ssb / dfb, ratio(ssw, dfw));
|
||
let f = msw.as_f64().map_or(f64::NAN, |w| msb / w);
|
||
g.txt(row + 2, 0, "Between Groups");
|
||
g.num(row + 2, 1, ssb);
|
||
g.num(row + 2, 2, dfb);
|
||
g.num(row + 2, 3, msb);
|
||
g.put(row + 2, 4, n_or_div0(f));
|
||
g.put(row + 2, 5, if f.is_finite() && dfw > 0.0 { callf("F.DIST.RT", &[f, dfb, dfw]) } else { div0() });
|
||
g.put(row + 2, 6, if dfw > 0.0 { callf("F.INV.RT", &[alpha, dfb, dfw]) } else { div0() });
|
||
g.txt(row + 3, 0, "Within Groups");
|
||
g.num(row + 3, 1, ssw);
|
||
g.num(row + 3, 2, dfw);
|
||
g.put(row + 3, 3, msw);
|
||
g.txt(row + 5, 0, "Total");
|
||
g.num(row + 5, 1, ssb + ssw);
|
||
g.num(row + 5, 2, n - 1.0);
|
||
Ok(g)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- two-factor anova and z-test
|
||
|
||
/// Count, Sum, Average and Variance of a group, as the Anova summaries list them.
|
||
fn group_summary(x: &[f64]) -> [Value; 4] {
|
||
[Value::number(x.len() as f64), Value::number(x.iter().sum()), n_or_div0(mean(x)), n_or_div0(var_s(x))]
|
||
}
|
||
|
||
const ANOVA_HEAD: [&str; 7] = ["Source of Variation", "SS", "df", "MS", "F", "P-value", "F crit"];
|
||
|
||
/// One line of an Anova table. A factor line, tested against the error term `err` (its MS and
|
||
/// df), adds F, the P-value and the critical F.
|
||
fn anova_line(g: &mut Grid, row: u32, label: &str, ss: f64, df: f64, err: Option<(f64, f64)>, alpha: f64) {
|
||
g.txt(row, 0, label);
|
||
g.num(row, 1, ss);
|
||
g.num(row, 2, df);
|
||
g.put(row, 3, ratio(ss, df));
|
||
if let Some((mse, dfe)) = err {
|
||
let f = ss / df / mse;
|
||
g.put(row, 4, n_or_div0(f));
|
||
g.put(row, 5, if f.is_finite() && dfe > 0.0 { callf("F.DIST.RT", &[f, df, dfe]) } else { div0() });
|
||
g.put(row, 6, if df > 0.0 && dfe > 0.0 { callf("F.INV.RT", &[alpha, df, dfe]) } else { div0() });
|
||
}
|
||
}
|
||
|
||
/// The cells of a two-factor input range, row by row.
|
||
fn input_rows(s: &Session, p: &Json) -> Result<Vec<Vec<Value>>> {
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
Ok(variables(s, si, r, true, false, &|_| String::new())?.into_iter().map(|v| v.vals).collect())
|
||
}
|
||
|
||
/// The values of column `j` of a table of numbers.
|
||
fn column(data: &[Vec<f64>], j: usize) -> Vec<f64> {
|
||
data.iter().filter_map(|r| r.get(j).copied()).collect()
|
||
}
|
||
|
||
/// Anova: Two-Factor Without Replication: one observation for each pair of levels (rows ×
|
||
/// columns), the error term being the residual after both main effects.
|
||
fn anova_two_factor_without(s: &Session, p: &Json) -> Result<Grid> {
|
||
let alpha = alpha(p)?;
|
||
let labels = flag(p, "labels");
|
||
let mut rows = input_rows(s, p)?;
|
||
let head = if labels && !rows.is_empty() { Some(rows.remove(0)) } else { None };
|
||
let mut row_names = Vec::with_capacity(rows.len());
|
||
let mut data: Vec<Vec<f64>> = Vec::with_capacity(rows.len());
|
||
for (i, mut row) in rows.into_iter().enumerate() {
|
||
let name = if labels && !row.is_empty() { row.remove(0).display() } else { format!("Row {}", i + 1) };
|
||
row_names.push(name);
|
||
data.push(numbers(&row, false)?);
|
||
}
|
||
let (a, b) = (data.len(), data.first().map_or(0, Vec::len));
|
||
if a < 2 || b < 2 {
|
||
return Err(EngineError::Other("Anova needs at least two rows and two columns of numbers.".into()));
|
||
}
|
||
let col_names: Vec<String> = match head {
|
||
Some(h) => h.iter().skip(1).map(Value::display).collect(),
|
||
None => (1..=b).map(|j| format!("Column {j}")).collect(),
|
||
};
|
||
let cols: Vec<Vec<f64>> = (0..b).map(|j| column(&data, j)).collect();
|
||
let all: Vec<f64> = data.iter().flatten().copied().collect();
|
||
let grand = mean(&all);
|
||
let row_means: Vec<f64> = data.iter().map(|r| mean(r)).collect();
|
||
let col_means: Vec<f64> = cols.iter().map(|c| mean(c)).collect();
|
||
let ss_rows: f64 = row_means.iter().map(|m| b as f64 * (m - grand).powi(2)).sum();
|
||
let ss_cols: f64 = col_means.iter().map(|m| a as f64 * (m - grand).powi(2)).sum();
|
||
// The residual of each observation after removing both main effects.
|
||
let ss_err: f64 =
|
||
data.iter().zip(&row_means).flat_map(|(r, rm)| r.iter().zip(&col_means).map(move |(x, cm)| (x - rm - cm + grand).powi(2))).sum();
|
||
let (df_rows, df_cols) = ((a - 1) as f64, (b - 1) as f64);
|
||
let df_err = df_rows * df_cols;
|
||
let mut g = Grid::default();
|
||
g.txt(0, 0, "Anova: Two-Factor Without Replication");
|
||
g.labels(2, 0, &["SUMMARY", "Count", "Sum", "Average", "Variance"]);
|
||
let summary = |g: &mut Grid, row: u32, name: &str, x: &[f64]| {
|
||
g.txt(row, 0, name);
|
||
for (k, v) in group_summary(x).into_iter().enumerate() {
|
||
g.put(row, 1 + idx(k), v);
|
||
}
|
||
};
|
||
for (i, (name, x)) in row_names.iter().zip(&data).enumerate() {
|
||
summary(&mut g, 3 + idx(i), name, x);
|
||
}
|
||
let first_col = 3 + idx(a) + 1;
|
||
for (j, (name, x)) in col_names.iter().zip(&cols).enumerate() {
|
||
summary(&mut g, first_col + idx(j), name, x);
|
||
}
|
||
let row = first_col + idx(b) + 2;
|
||
g.txt(row, 0, "ANOVA");
|
||
g.labels(row + 1, 0, &ANOVA_HEAD);
|
||
let err = Some((ss_err / df_err, df_err));
|
||
anova_line(&mut g, row + 2, "Rows", ss_rows, df_rows, err, alpha);
|
||
anova_line(&mut g, row + 3, "Columns", ss_cols, df_cols, err, alpha);
|
||
anova_line(&mut g, row + 4, "Error", ss_err, df_err, None, alpha);
|
||
g.txt(row + 6, 0, "Total");
|
||
g.num(row + 6, 1, ss(&all));
|
||
g.num(row + 6, 2, (all.len() - 1) as f64);
|
||
Ok(g)
|
||
}
|
||
|
||
/// The number of rows in each sample of a two-factor Anova with replication.
|
||
fn rows_per_sample(p: &Json, data_rows: usize) -> Result<usize> {
|
||
let per = num(p, "rowsPerSample").ok_or_else(|| bad(CMD, "missing `rowsPerSample` (the number of rows in each sample)"))?;
|
||
if per < 2.0 || per.fract() != 0.0 || per > data_rows as f64 || !data_rows.is_multiple_of(per as usize) {
|
||
return Err(EngineError::Other("Rows per sample must be a whole number of at least 2 that divides the number of data rows.".into()));
|
||
}
|
||
Ok(per as usize)
|
||
}
|
||
|
||
/// Anova: Two-Factor With Replication. The input range holds column labels in its first row and
|
||
/// sample labels in its first column; each sample (a level of the row factor) spans
|
||
/// `rowsPerSample` consecutive rows, so every cell of the design holds that many observations.
|
||
fn anova_two_factor_with(s: &Session, p: &Json) -> Result<Grid> {
|
||
let alpha = alpha(p)?;
|
||
let mut rows = input_rows(s, p)?;
|
||
if rows.len() < 3 {
|
||
return Err(EngineError::Other("The input range must hold a row of labels and at least two samples.".into()));
|
||
}
|
||
let head = rows.remove(0);
|
||
let col_names: Vec<String> = head.iter().skip(1).map(Value::display).collect();
|
||
let per = rows_per_sample(p, rows.len())?;
|
||
let mut sample_names = Vec::new();
|
||
let mut data: Vec<Vec<f64>> = Vec::with_capacity(rows.len());
|
||
for (i, mut row) in rows.into_iter().enumerate() {
|
||
let name = if row.is_empty() { Value::Empty } else { row.remove(0) };
|
||
if i.is_multiple_of(per) {
|
||
sample_names.push(name.display());
|
||
}
|
||
data.push(numbers(&row, false)?);
|
||
}
|
||
let (a, b) = (sample_names.len(), col_names.len());
|
||
if a < 2 || b < 2 {
|
||
return Err(EngineError::Other("Anova needs at least two samples and two columns of numbers.".into()));
|
||
}
|
||
// cells[i][j]: the observations of sample i in column j.
|
||
let cells: Vec<Vec<Vec<f64>>> = data.chunks(per).map(|block| (0..b).map(|j| column(block, j)).collect()).collect();
|
||
let all: Vec<f64> = data.iter().flatten().copied().collect();
|
||
let grand = mean(&all);
|
||
let sample_all: Vec<Vec<f64>> = data.chunks(per).map(|block| block.iter().flatten().copied().collect()).collect();
|
||
let cols: Vec<Vec<f64>> = (0..b).map(|j| column(&data, j)).collect();
|
||
let (sm, cm): (Vec<f64>, Vec<f64>) = (sample_all.iter().map(|x| mean(x)).collect(), cols.iter().map(|x| mean(x)).collect());
|
||
let n = per as f64;
|
||
let ss_sample: f64 = sm.iter().map(|m| n * b as f64 * (m - grand).powi(2)).sum();
|
||
let ss_cols: f64 = cm.iter().map(|m| n * a as f64 * (m - grand).powi(2)).sum();
|
||
let mut ss_inter = 0.0;
|
||
let mut ss_within = 0.0;
|
||
for (i, row) in cells.iter().enumerate() {
|
||
for (j, cell) in row.iter().enumerate() {
|
||
let (si, cj) = (sm.get(i).copied().unwrap_or(f64::NAN), cm.get(j).copied().unwrap_or(f64::NAN));
|
||
ss_inter += n * (mean(cell) - si - cj + grand).powi(2);
|
||
ss_within += ss(cell);
|
||
}
|
||
}
|
||
let (df_sample, df_cols) = ((a - 1) as f64, (b - 1) as f64);
|
||
let df_within = (a * b * (per - 1)) as f64;
|
||
let mut g = Grid::default();
|
||
g.txt(0, 0, "Anova: Two-Factor With Replication");
|
||
g.txt(2, 0, "SUMMARY");
|
||
for (j, name) in col_names.iter().enumerate() {
|
||
g.txt(2, 1 + idx(j), name);
|
||
}
|
||
g.txt(2, 1 + idx(b), "Total");
|
||
let stats = ["Count", "Sum", "Average", "Variance"];
|
||
let mut row = 3u32;
|
||
let block = |g: &mut Grid, row: &mut u32, name: &str, groups: Vec<&[f64]>| {
|
||
g.txt(*row, 0, name);
|
||
let sums: Vec<[Value; 4]> = groups.iter().map(|x| group_summary(x)).collect();
|
||
for (k, label) in stats.iter().enumerate() {
|
||
let r = *row + 1 + idx(k);
|
||
g.txt(r, 0, label);
|
||
for (j, sm) in sums.iter().enumerate() {
|
||
g.put(r, 1 + idx(j), sm.get(k).cloned().unwrap_or(Value::Empty));
|
||
}
|
||
}
|
||
*row += 6;
|
||
};
|
||
for ((name, row_cells), whole) in sample_names.iter().zip(&cells).zip(&sample_all) {
|
||
let mut groups: Vec<&[f64]> = row_cells.iter().map(Vec::as_slice).collect();
|
||
groups.push(whole);
|
||
block(&mut g, &mut row, name, groups);
|
||
}
|
||
block(&mut g, &mut row, "Total", cols.iter().map(Vec::as_slice).collect());
|
||
row += 1;
|
||
g.txt(row, 0, "ANOVA");
|
||
g.labels(row + 1, 0, &ANOVA_HEAD);
|
||
let err = Some((ss_within / df_within, df_within));
|
||
anova_line(&mut g, row + 2, "Sample", ss_sample, df_sample, err, alpha);
|
||
anova_line(&mut g, row + 3, "Columns", ss_cols, df_cols, err, alpha);
|
||
anova_line(&mut g, row + 4, "Interaction", ss_inter, df_sample * df_cols, err, alpha);
|
||
anova_line(&mut g, row + 5, "Within", ss_within, df_within, None, alpha);
|
||
g.txt(row + 7, 0, "Total");
|
||
g.num(row + 7, 1, ss(&all));
|
||
g.num(row + 7, 2, (all.len() - 1) as f64);
|
||
Ok(g)
|
||
}
|
||
|
||
/// A known variance parameter: a positive number.
|
||
fn known_variance(p: &Json, key: &str) -> Result<f64> {
|
||
match num(p, key) {
|
||
Some(v) if v > 0.0 => Ok(v),
|
||
_ => Err(EngineError::Other(format!("Enter the known variance `{key}` (a positive number)."))),
|
||
}
|
||
}
|
||
|
||
/// z-Test: Two Sample for Means, for two samples whose population variances are known.
|
||
fn z_test(s: &Session, p: &Json) -> Result<Grid> {
|
||
let alpha = alpha(p)?;
|
||
let hmd = match opt_num(p, "hypothesizedDifference")? {
|
||
Some(d) => d,
|
||
None => num_or(p, "difference", 0.0)?,
|
||
};
|
||
let (v1, v2) = (known_variance(p, "variance1")?, known_variance(p, "variance2")?);
|
||
let a = single(s, p, "range1", "Variable 1")?;
|
||
let b = single(s, p, "range2", "Variable 2")?;
|
||
let (x, y) = (numbers(&a.vals, true)?, numbers(&b.vals, true)?);
|
||
if x.is_empty() || y.is_empty() {
|
||
return Err(EngineError::Other("Each input range needs at least one number.".into()));
|
||
}
|
||
let (n1, n2) = (x.len() as f64, y.len() as f64);
|
||
let (m1, m2) = (mean(&x), mean(&y));
|
||
let z = (m1 - m2 - hmd) / (v1 / n1 + v2 / n2).sqrt();
|
||
let mut g = Grid::default();
|
||
let row = test_head(&mut g, "z-Test: Two Sample for Means", &a.name, &b.name);
|
||
pair_row(&mut g, row, "Mean", Value::number(m1), Value::number(m2));
|
||
pair_row(&mut g, row + 1, "Known Variance", Value::number(v1), Value::number(v2));
|
||
pair_row(&mut g, row + 2, "Observations", Value::number(n1), Value::number(n2));
|
||
pair_row(&mut g, row + 3, "Hypothesized Mean Difference", Value::number(hmd), Value::Empty);
|
||
pair_row(&mut g, row + 4, "z", n_or_div0(z), Value::Empty);
|
||
// The tail beyond |z|: 1 − Φ(|z|).
|
||
let tail = call("NORM.S.DIST", &[Value::Number(z.abs()), Value::Bool(true)]).as_f64().map(|c| 1.0 - c);
|
||
let lines: [(&str, Value); 4] = [
|
||
("P(Z<=z) one-tail", tail.map_or_else(div0, Value::number)),
|
||
("z Critical one-tail", callf("NORM.S.INV", &[1.0 - alpha])),
|
||
("P(Z<=z) two-tail", tail.map_or_else(div0, |t| Value::number(2.0 * t))),
|
||
("z Critical two-tail", callf("NORM.S.INV", &[1.0 - alpha / 2.0])),
|
||
];
|
||
for (k, (label, v)) in lines.into_iter().enumerate() {
|
||
pair_row(&mut g, row + 5 + idx(k), label, v, Value::Empty);
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- fourier analysis
|
||
|
||
/// Radix-2 fast Fourier transform (iterative Cooley–Tukey) in place, for a power-of-two number of
|
||
/// complex values (re, im). The forward transform is unscaled, X[k] = Σ x[n]·e^(−2πikn/N); the
|
||
/// inverse uses e^(+2πikn/N) and divides by N, so it returns the values the forward one took.
|
||
fn fft(a: &mut [(f64, f64)], inverse: bool) {
|
||
let n = a.len();
|
||
if n < 2 || !n.is_power_of_two() {
|
||
return;
|
||
}
|
||
let bits = n.trailing_zeros();
|
||
for i in 0..n {
|
||
let j = i.reverse_bits() >> (usize::BITS - bits);
|
||
if i < j {
|
||
a.swap(i, j);
|
||
}
|
||
}
|
||
let sign = if inverse { 1.0 } else { -1.0 };
|
||
let mut len = 2;
|
||
while len <= n {
|
||
let ang = sign * 2.0 * std::f64::consts::PI / len as f64;
|
||
for start in (0..n).step_by(len) {
|
||
for k in 0..len / 2 {
|
||
// Each twiddle factor is computed directly, not by repeated multiplication.
|
||
let (wr, wi) = ((ang * k as f64).cos(), (ang * k as f64).sin());
|
||
let (i, j) = (start + k, start + k + len / 2);
|
||
let (Some(&(ur, ui)), Some(&(xr, xi))) = (a.get(i), a.get(j)) else { continue };
|
||
let (vr, vi) = (xr * wr - xi * wi, xr * wi + xi * wr);
|
||
if let Some(u) = a.get_mut(i) {
|
||
*u = (ur + vr, ui + vi);
|
||
}
|
||
if let Some(x) = a.get_mut(j) {
|
||
*x = (ur - vr, ui - vi);
|
||
}
|
||
}
|
||
}
|
||
len <<= 1;
|
||
}
|
||
if inverse {
|
||
let nf = n as f64;
|
||
a.iter_mut().for_each(|v| *v = (v.0 / nf, v.1 / nf));
|
||
}
|
||
}
|
||
|
||
/// Fourier Analysis: the discrete Fourier transform (or with `inverse` its inverse) of one row or
|
||
/// column of 2^k values, read and written as complex numbers the way IMREAL, IMAGINARY and
|
||
/// COMPLEX read and write them.
|
||
fn fourier_analysis(s: &Session, p: &Json) -> Result<Grid> {
|
||
let v = single(s, p, "range", "")?;
|
||
let n = v.vals.len();
|
||
if n == 0 || n > MAX_FOURIER || !n.is_power_of_two() {
|
||
return Err(EngineError::Other(format!("Fourier Analysis needs 1, 2, 4, 8 … {MAX_FOURIER} input values (a power of two).")));
|
||
}
|
||
let mut a = Vec::with_capacity(n);
|
||
for x in &v.vals {
|
||
let part = |f: &str| call(f, std::slice::from_ref(x)).as_f64();
|
||
match (x, part("IMREAL"), part("IMAGINARY")) {
|
||
(Value::Empty, ..) => return Err(non_numeric()),
|
||
(_, Some(re), Some(im)) => a.push((re, im)),
|
||
_ => return Err(EngineError::Other("Input range must hold numbers or complex numbers such as 3+4i.".into())),
|
||
}
|
||
}
|
||
let total: f64 = a.iter().map(|(re, im)| re.abs() + im.abs()).sum();
|
||
let inverse = flag(p, "inverse");
|
||
fft(&mut a, inverse);
|
||
// Parts that can only be rounding noise of the transform are written as 0. Each of the
|
||
// log2(N) butterfly stages adds at most a few units of f64::EPSILON times the magnitude of
|
||
// the values it combines, so the computed transform differs from the exact one by at most
|
||
// about log2(N) · ε · Σ|x_n| times a small constant (Higham, Accuracy and Stability of
|
||
// Numerical Algorithms, 2nd ed., chapter 24, The Fast Fourier Transform). The bound comes
|
||
// from the input and N, never from the largest output, so a small part beside a large one
|
||
// (an imaginary part of 1 beside a real part of 1e14, a harmonic of 4e-7 beside a DC term of
|
||
// 8e6) is kept. With N = 1 the transform is the identity and log2(N) = 0: nothing is zeroed.
|
||
// The inverse divides its outputs by N, so its noise is divided by N as well. An input sum
|
||
// that overflows gives no finite bound; then nothing is cleaned up.
|
||
let mut noise = 8.0 * f64::EPSILON * f64::from(n.trailing_zeros()) * total;
|
||
if inverse {
|
||
noise /= n as f64;
|
||
}
|
||
let noise = if noise.is_finite() { noise } else { 0.0 };
|
||
let clean = |x: f64| if x.abs() <= noise { 0.0 } else { x };
|
||
let mut g = Grid::default();
|
||
let first = u32::from(flag(p, "labels"));
|
||
if first == 1 {
|
||
g.txt(0, 0, &v.name);
|
||
}
|
||
for (i, (re, im)) in a.into_iter().enumerate() {
|
||
g.put(first + idx(i), 0, callf("COMPLEX", &[clean(re), clean(im)]));
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- random numbers and sampling
|
||
|
||
/// A seeded generator: xoshiro256** (Blackman and Vigna, a public-domain algorithm) with its
|
||
/// state filled by SplitMix64. The same seed gives the same numbers in GridCraft; they are not
|
||
/// Excel's numbers, whose generator is not documented.
|
||
struct Rng([u64; 4]);
|
||
|
||
impl Rng {
|
||
fn new(seed: u64) -> Rng {
|
||
let mut x = seed;
|
||
let mut next = || {
|
||
x = x.wrapping_add(0x9E37_79B9_7F4A_7C15);
|
||
let mut z = x;
|
||
z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
|
||
z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
|
||
z ^ (z >> 31)
|
||
};
|
||
Rng([next(), next(), next(), next()])
|
||
}
|
||
|
||
fn next_u64(&mut self) -> u64 {
|
||
let [s0, s1, s2, s3] = self.0;
|
||
let out = s1.wrapping_mul(5).rotate_left(7).wrapping_mul(9);
|
||
let (s2, s3) = (s2 ^ s0, s3 ^ s1);
|
||
self.0 = [s0 ^ s3, s1 ^ s2, s2 ^ (s1 << 17), s3.rotate_left(45)];
|
||
out
|
||
}
|
||
|
||
/// Uniform in [0, 1), with 53 random bits.
|
||
fn uniform(&mut self) -> f64 {
|
||
(self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
|
||
}
|
||
|
||
/// Uniform whole number in [0, n), without modulo bias.
|
||
fn below(&mut self, n: u64) -> u64 {
|
||
if n <= 1 {
|
||
return 0;
|
||
}
|
||
let zone = u64::MAX - u64::MAX % n;
|
||
loop {
|
||
let v = self.next_u64();
|
||
if v < zone {
|
||
return v % n;
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Standard normal by the Box–Muller transform.
|
||
fn normal(&mut self) -> f64 {
|
||
let (u1, u2) = (1.0 - self.uniform(), self.uniform());
|
||
(-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
|
||
}
|
||
}
|
||
|
||
/// The `seed` parameter, or a seed taken from the clock when there is none.
|
||
fn seed_of(p: &Json) -> Result<u64> {
|
||
let given = p.get("seed").is_some_and(|v| match v {
|
||
Json::Null => false,
|
||
Json::String(t) => !t.trim().is_empty(),
|
||
_ => true,
|
||
});
|
||
if !given {
|
||
let ms = gridcraft_core::date::unix_now_ms().unwrap_or(0.0);
|
||
return Ok((ms * 1000.0) as u64 & MAX_SEED);
|
||
}
|
||
match num(p, "seed") {
|
||
Some(x) if x >= 0.0 && x.fract() == 0.0 && x <= MAX_SEED as f64 => Ok(x as u64),
|
||
_ => Err(EngineError::Other("The random seed must be a whole number from 0.".into())),
|
||
}
|
||
}
|
||
|
||
/// A number parameter: `None` when absent, null or blank; an error when present but not a
|
||
/// finite number (a typed "0,3" or "abc" must not quietly become the default).
|
||
fn opt_num(p: &Json, key: &str) -> Result<Option<f64>> {
|
||
match p.get(key) {
|
||
None | Some(Json::Null) => Ok(None),
|
||
Some(Json::String(t)) if t.trim().is_empty() => Ok(None),
|
||
Some(_) => num(p, key).map(Some).ok_or_else(|| EngineError::Other(format!("`{key}` must be a number."))),
|
||
}
|
||
}
|
||
|
||
/// A number parameter, `default` when absent or blank.
|
||
fn num_or(p: &Json, key: &str, default: f64) -> Result<f64> {
|
||
Ok(opt_num(p, key)?.unwrap_or(default))
|
||
}
|
||
|
||
/// A whole-number parameter of at least `min`, `default` when absent.
|
||
fn whole(p: &Json, key: &str, default: Option<f64>, min: f64) -> Result<f64> {
|
||
match opt_num(p, key)?.or(default) {
|
||
Some(x) if x >= min && x.fract() == 0.0 => Ok(x),
|
||
Some(_) => Err(EngineError::Other(format!("`{key}` must be a whole number of at least {min}."))),
|
||
None => Err(bad(CMD, format!("missing `{key}`"))),
|
||
}
|
||
}
|
||
|
||
/// A probability parameter in [0, 1].
|
||
fn probability(p: &Json) -> Result<f64> {
|
||
match num_or(p, "p", 0.5)? {
|
||
x if (0.0..=1.0).contains(&x) => Ok(x),
|
||
_ => Err(EngineError::Other("The p value must be between 0 and 1.".into())),
|
||
}
|
||
}
|
||
|
||
/// The probabilities of a distribution on 0, 1, 2 … `last`, as (first outcome, probabilities of
|
||
/// first, first + 1, …), from its probability `at_mode` at `mode`, its most likely outcome, and
|
||
/// `ratio`(k) = pmf(k + 1) / pmf(k). Each term is one multiplication of its neighbour towards the
|
||
/// mode, so no worksheet function runs per outcome, and the terms only shrink away from the mode,
|
||
/// so nothing underflows before the terms that matter. The left tail ends where its terms vanish,
|
||
/// the right tail once the cumulative probability reaches 1 within f64 resolution (every later
|
||
/// draw would map to the last outcome). The table is as wide as the distribution's effective
|
||
/// support, not as `last`: Binomial(1,000,000, 0.5) takes about 23,000 entries, not 1,000,001.
|
||
fn pmf_table(last: f64, mode: f64, at_mode: f64, ratio: impl Fn(f64) -> f64) -> Result<(f64, Vec<f64>)> {
|
||
if !at_mode.is_finite() || at_mode <= 0.0 || !(0.0..=last).contains(&mode) || mode.fract() != 0.0 {
|
||
return Err(EngineError::Other("The distribution could not be tabulated.".into()));
|
||
}
|
||
let mut left = Vec::new();
|
||
let (mut k, mut v) = (mode, at_mode);
|
||
while k > 0.0 {
|
||
v /= ratio(k - 1.0);
|
||
if v.is_nan() || v < f64::MIN_POSITIVE {
|
||
break;
|
||
}
|
||
left.push(v);
|
||
k -= 1.0;
|
||
}
|
||
let first = k;
|
||
let mut pmf: Vec<f64> = left.into_iter().rev().collect();
|
||
let mut acc = pmf.iter().sum::<f64>() + at_mode;
|
||
pmf.push(at_mode);
|
||
let (mut k, mut v) = (mode, at_mode);
|
||
while k < last && acc < 1.0 - 1e-15 {
|
||
v *= ratio(k);
|
||
if v.is_nan() || v < f64::MIN_POSITIVE {
|
||
break;
|
||
}
|
||
pmf.push(v);
|
||
acc += v;
|
||
k += 1.0;
|
||
}
|
||
Ok((first, pmf))
|
||
}
|
||
|
||
/// A distribution on whole numbers tabulated for inversion sampling: the cumulative
|
||
/// probabilities of `first`, `first + 1`, … up to the last outcome that matters.
|
||
struct Table {
|
||
first: f64,
|
||
cdf: Vec<f64>,
|
||
}
|
||
|
||
impl Table {
|
||
/// The distribution on 0, 1, 2 … `last` that `pmf_table` describes.
|
||
fn from_mode(last: f64, mode: f64, at_mode: f64, ratio: impl Fn(f64) -> f64) -> Result<Table> {
|
||
let (first, pmf) = pmf_table(last, mode, at_mode, ratio)?;
|
||
let mut acc = 0.0;
|
||
let cdf = pmf
|
||
.into_iter()
|
||
.map(|p| {
|
||
acc += p;
|
||
acc
|
||
})
|
||
.collect();
|
||
Ok(Table { first, cdf })
|
||
}
|
||
|
||
/// All the probability on one outcome.
|
||
fn point(value: f64) -> Table {
|
||
Table { first: value, cdf: vec![1.0] }
|
||
}
|
||
|
||
/// A draw: the first outcome whose cumulative probability exceeds a uniform number (inversion).
|
||
fn draw(&self, r: &mut Rng) -> f64 {
|
||
let u = r.uniform() * self.cdf.last().copied().unwrap_or(1.0);
|
||
self.first + self.cdf.partition_point(|c| *c <= u).min(self.cdf.len().saturating_sub(1)) as f64
|
||
}
|
||
}
|
||
|
||
/// The numbers of a patterned sequence: `from` to `to` in steps, each number repeated, the whole
|
||
/// sequence repeated.
|
||
fn pattern(p: &Json) -> Result<Vec<f64>> {
|
||
let (from, to, step) = (num_or(p, "from", 1.0)?, num_or(p, "to", 1.0)?, num_or(p, "step", 1.0)?);
|
||
let (each, times) = (whole(p, "repeatEach", Some(1.0), 1.0)?, whole(p, "repeatSequence", Some(1.0), 1.0)?);
|
||
if step <= 0.0 || to < from {
|
||
return Err(EngineError::Other("The pattern must run from a smaller number to a larger one in positive steps.".into()));
|
||
}
|
||
let steps = ((to - from) / step + 1e-9).floor() + 1.0;
|
||
if steps * each * times > MAX_RANDOM as f64 {
|
||
return Err(EngineError::Other(format!("Random Number Generation writes at most {MAX_RANDOM} numbers.")));
|
||
}
|
||
let base: Vec<f64> = (0..steps as usize).map(|k| from + step * k as f64).collect();
|
||
let one: Vec<f64> = base.iter().flat_map(|v| std::iter::repeat_n(*v, each as usize)).collect();
|
||
Ok(one.repeat(times as usize))
|
||
}
|
||
|
||
/// Random Number Generation: `count` numbers from a distribution in each of `variables`
|
||
/// columns, drawn row by row.
|
||
fn random_numbers(s: &Session, p: &Json, seed: u64) -> Result<Grid> {
|
||
let vars = whole(p, "variables", Some(1.0), 1.0)?;
|
||
let dist = text(p, &["distribution"]).map(norm).unwrap_or_else(|| "uniform".into());
|
||
let cells = |rows: f64| {
|
||
if rows * vars > MAX_RANDOM as f64 {
|
||
Err(EngineError::Other(format!("Random Number Generation writes at most {MAX_RANDOM} numbers.")))
|
||
} else {
|
||
Ok(rows as usize)
|
||
}
|
||
};
|
||
let mut g = Grid::default();
|
||
let cols = vars as usize;
|
||
if dist == "patterned" {
|
||
let seq = pattern(p)?;
|
||
cells(seq.len() as f64)?;
|
||
for (i, x) in seq.iter().enumerate() {
|
||
for c in 0..cols {
|
||
g.num(idx(i), idx(c), *x);
|
||
}
|
||
}
|
||
return Ok(g);
|
||
}
|
||
let count = cells(whole(p, "count", None, 1.0)?)?;
|
||
let mut draw: Box<dyn FnMut(&mut Rng) -> f64> = match dist.as_str() {
|
||
"uniform" => {
|
||
let (a, b) = (num_or(p, "from", 0.0)?, num_or(p, "to", 1.0)?);
|
||
Box::new(move |r| a + (b - a) * r.uniform())
|
||
}
|
||
"normal" => {
|
||
let (m, sd) = (num_or(p, "mean", 0.0)?, num_or(p, "standardDeviation", 1.0)?);
|
||
if sd < 0.0 {
|
||
return Err(EngineError::Other("The standard deviation must not be negative.".into()));
|
||
}
|
||
Box::new(move |r| m + sd * r.normal())
|
||
}
|
||
"bernoulli" => {
|
||
let pr = probability(p)?;
|
||
Box::new(move |r| if r.uniform() < pr { 1.0 } else { 0.0 })
|
||
}
|
||
"binomial" => {
|
||
let (pr, n) = (probability(p)?, whole(p, "trials", Some(1.0), 1.0)?);
|
||
if n > MAX_OUTCOMES {
|
||
return Err(EngineError::Other(format!("The number of trials must be at most {MAX_OUTCOMES}.")));
|
||
}
|
||
// p = 0 and p = 1 are exact: every draw is 0 or n. Otherwise one BINOM.DIST value at
|
||
// the mode ⌊(n + 1)p⌋ anchors pmf(k + 1) = pmf(k) · (n − k)/(k + 1) · p/(1 − p).
|
||
let table = if pr == 0.0 {
|
||
Table::point(0.0)
|
||
} else if pr == 1.0 {
|
||
Table::point(n)
|
||
} else {
|
||
let mode = ((n + 1.0) * pr).floor().min(n);
|
||
let at_mode = f64_of(&call("BINOM.DIST", &[Value::Number(mode), Value::Number(n), Value::Number(pr), Value::Bool(false)]));
|
||
let odds = pr / (1.0 - pr);
|
||
Table::from_mode(n, mode, at_mode, |k| (n - k) / (k + 1.0) * odds)?
|
||
};
|
||
Box::new(move |r| table.draw(r))
|
||
}
|
||
"poisson" => {
|
||
let lambda = num_or(p, "lambda", 1.0)?;
|
||
if lambda <= 0.0 || lambda > 100_000.0 {
|
||
return Err(EngineError::Other("Lambda must be greater than 0 and at most 100000.".into()));
|
||
}
|
||
// One POISSON.DIST value at the mode ⌊λ⌋ anchors pmf(k + 1) = pmf(k) · λ/(k + 1).
|
||
let last = (lambda + 40.0 * lambda.sqrt() + 100.0).min(MAX_OUTCOMES);
|
||
let mode = lambda.floor();
|
||
let at_mode = f64_of(&call("POISSON.DIST", &[Value::Number(mode), Value::Number(lambda), Value::Bool(false)]));
|
||
let table = Table::from_mode(last, mode, at_mode, |k| lambda / (k + 1.0))?;
|
||
Box::new(move |r| table.draw(r))
|
||
}
|
||
"discrete" => {
|
||
let (si, r) = input_range(s, p, &["valueRange", "values"])?;
|
||
if r.width() != 2 {
|
||
return Err(EngineError::Other("The value and probability input range must have two columns.".into()));
|
||
}
|
||
let v = variables(s, si, r, false, false, &|_| String::new())?;
|
||
let (values, probs) = match v.as_slice() {
|
||
[a, b] => (numbers(&a.vals, false)?, numbers(&b.vals, false)?),
|
||
_ => return Err(non_numeric()),
|
||
};
|
||
if probs.iter().any(|x| !(0.0..=1.0).contains(x)) || (probs.iter().sum::<f64>() - 1.0).abs() > 1e-9 {
|
||
return Err(EngineError::Other("The probabilities must be between 0 and 1 and add up to 1.".into()));
|
||
}
|
||
let cdf: Vec<f64> = probs
|
||
.iter()
|
||
.scan(0.0, |acc, x| {
|
||
*acc += x;
|
||
Some(*acc)
|
||
})
|
||
.collect();
|
||
let table = Table { first: 0.0, cdf };
|
||
Box::new(move |r| values.get(table.draw(r) as usize).copied().unwrap_or(f64::NAN))
|
||
}
|
||
other => return Err(bad(CMD, format!("unknown distribution `{other}`"))),
|
||
};
|
||
let mut rng = Rng::new(seed);
|
||
for i in 0..count {
|
||
for c in 0..cols {
|
||
g.num(idx(i), idx(c), draw(&mut rng));
|
||
}
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
/// Whether Random Number Generation writes a patterned sequence (which uses no seed).
|
||
fn patterned(p: &Json) -> bool {
|
||
text(p, &["distribution"]).map(norm).is_some_and(|d| d == "patterned")
|
||
}
|
||
|
||
fn random_sampling(p: &Json) -> bool {
|
||
text(p, &["method"]).is_some_and(|m| m.eq_ignore_ascii_case("random"))
|
||
}
|
||
|
||
/// Sampling: every `period`-th value of the input, or `samples` values drawn at random with
|
||
/// replacement, down one column. The input is read column by column (row by row for a single
|
||
/// row); blanks are skipped.
|
||
fn sampling(s: &Session, p: &Json, seed: Option<u64>) -> Result<Grid> {
|
||
let (si, r) = input_range(s, p, &["range", "input"])?;
|
||
let by_rows = r.height() == 1 && r.width() > 1;
|
||
let vars = variables(s, si, r, by_rows, flag(p, "labels"), &|_| String::new())?;
|
||
let mut x = Vec::new();
|
||
for v in &vars {
|
||
x.extend(numbers(&v.vals, true)?);
|
||
}
|
||
if x.is_empty() {
|
||
return Err(EngineError::Other("The input range holds no numbers.".into()));
|
||
}
|
||
let picked: Vec<f64> = match seed {
|
||
Some(seed) => {
|
||
let k = whole(p, "samples", None, 1.0)?;
|
||
if k > MAX_RANDOM as f64 {
|
||
return Err(EngineError::Other(format!("Sampling takes at most {MAX_RANDOM} samples.")));
|
||
}
|
||
let mut rng = Rng::new(seed);
|
||
(0..k as usize).filter_map(|_| x.get(rng.below(x.len() as u64) as usize).copied()).collect()
|
||
}
|
||
None => {
|
||
let period = whole(p, "period", None, 1.0)?;
|
||
if period > x.len() as f64 {
|
||
return Err(EngineError::Other("The period must be at most the number of input values.".into()));
|
||
}
|
||
x.iter().skip(period as usize - 1).step_by(period as usize).copied().collect()
|
||
}
|
||
};
|
||
let mut g = Grid::default();
|
||
for (i, v) in picked.iter().enumerate() {
|
||
g.num(idx(i), 0, *v);
|
||
}
|
||
Ok(g)
|
||
}
|
||
|
||
// ---------------------------------------------------------------- output
|
||
|
||
/// Writes the report at `output` (refusing to overwrite data unless `overwrite`) or on a new
|
||
/// sheet.
|
||
fn write_report(s: &mut Session, p: &Json, tool: &str, g: Grid) -> Result<Json> {
|
||
let (rows, cols) = g.extent();
|
||
let (rows, cols) = (rows.max(1), cols.max(1));
|
||
if u64::from(rows) * u64::from(cols) > MAX_OUTPUT_CELLS {
|
||
return Err(bad(CMD, "the report is too large"));
|
||
}
|
||
let d = s.doc()?;
|
||
let target = match str_param(p, "outputTo").map(str::to_ascii_lowercase).as_deref() {
|
||
Some("newsheet") => None,
|
||
Some("range") => Some(text(p, &["output", "outputRange"]).ok_or_else(|| bad(CMD, "missing `output` (the output range)"))?),
|
||
_ => text(p, &["output", "outputRange"]),
|
||
};
|
||
let (si, at, new_sheet) = match target {
|
||
Some(t) => {
|
||
let (si, r) = parse_range(&d.wb, t, CMD)?;
|
||
(si, r.start, None)
|
||
}
|
||
None => {
|
||
if d.wb.protected_structure {
|
||
return Err(EngineError::Other("The workbook is protected and cannot be changed.".into()));
|
||
}
|
||
let name = match p.get("newSheet") {
|
||
Some(Json::String(n)) if !n.trim().is_empty() => n.trim().to_string(),
|
||
_ => d.wb.next_sheet_name(),
|
||
};
|
||
d.wb.check_sheet_name(&name, None).map_err(EngineError::Other)?;
|
||
(d.wb.sheets.len(), CellRef::new(0, 0), Some(name))
|
||
}
|
||
};
|
||
let end_row = at.row.checked_add(rows - 1).filter(|r| *r < MAX_ROWS);
|
||
let end_col = at.col.checked_add(cols - 1).filter(|c| *c < MAX_COLS);
|
||
let (Some(end_row), Some(end_col)) = (end_row, end_col) else {
|
||
return Err(EngineError::Other("The report doesn't fit on the sheet below and right of the output cell.".into()));
|
||
};
|
||
let rect = RangeRef::new(at, CellRef::new(end_row, end_col));
|
||
// Existing content inside the report's rectangle.
|
||
let mut occupied: Vec<CellRef> = Vec::new();
|
||
if new_sheet.is_none()
|
||
&& let Some(sh) = d.wb.sheet(si)
|
||
{
|
||
for row in rect.start.row..=rect.end.row {
|
||
for (c, cell) in sh.cells.row(row, rect.start.col, rect.end.col) {
|
||
if cell.has_content() {
|
||
occupied.push(CellRef::new(row, c));
|
||
}
|
||
}
|
||
}
|
||
if sh.spill.keys().any(|c| rect.contains(*c)) && occupied.is_empty() {
|
||
occupied.push(rect.start);
|
||
}
|
||
}
|
||
if !occupied.is_empty() && !flag(p, "overwrite") {
|
||
return Err(EngineError::Other(format!("The output range {} {OVERWRITE_ERROR}; pass overwrite: true to replace it.", rect.a1())));
|
||
}
|
||
let sheet_name = new_sheet.clone().or_else(|| d.wb.sheet(si).map(|sh| sh.name.clone())).unwrap_or_default();
|
||
let new_index = edit(s, |cx| {
|
||
let pct_fmt = NumFmt::new("0.00%");
|
||
if let Some(name) = &new_sheet {
|
||
cx.wb.sheets.push(std::sync::Arc::new(Sheet::new(name.clone())));
|
||
cx.structural = true;
|
||
}
|
||
let pct_styles: Vec<(CellRef, gridcraft_model::StyleId)> = {
|
||
let sh = cx.wb.sheet(si).ok_or(EngineError::NoDocument)?;
|
||
g.pct.iter().map(|(r, c)| CellRef::new(at.row + r, at.col + c)).map(|c| (c, sh.style_id(c))).collect()
|
||
};
|
||
let pct_styles: HashMap<CellRef, gridcraft_model::StyleId> =
|
||
pct_styles.into_iter().map(|(c, st)| (c, cx.wb.styles.derive(st, |s| s.num_fmt = pct_fmt.clone()))).collect();
|
||
let sh = cx.sheet_mut(si)?;
|
||
for c in &occupied {
|
||
let style = sh.style_id(*c);
|
||
sh.set_cell(*c, Cell { value: Value::Empty, formula: None, style, py: None });
|
||
}
|
||
let mut touched = occupied.clone();
|
||
for ((r, c), v) in &g.cells {
|
||
let cell = CellRef::new(at.row + r, at.col + c);
|
||
let style = pct_styles.get(&cell).copied().unwrap_or_else(|| sh.style_id(cell));
|
||
sh.set_cell(cell, Cell { value: v.clone(), formula: None, style, py: None });
|
||
touched.push(cell);
|
||
}
|
||
for c in touched {
|
||
cx.touch(si, c);
|
||
}
|
||
Ok(new_sheet.as_ref().map(|_| si))
|
||
})?;
|
||
if let Some(i) = new_index {
|
||
s.doc_mut()?.wb_switch_sheet(i);
|
||
}
|
||
Ok(json!({"tool": tool, "sheet": sheet_name, "range": rect.a1(), "rows": rows, "cols": cols}))
|
||
}
|
||
|
||
fn run(s: &mut Session, p: &Json) -> Result<Json> {
|
||
let Some(t) = text(p, &["tool"]) else {
|
||
s.ui_requests.push(crate::UiRequest::Dialog("analysisToolPak".into(), json!({})));
|
||
return Ok(json!({"tools": TOOLS.iter().map(|(id, label)| json!({"id": id, "label": label})).collect::<Vec<_>>()}));
|
||
};
|
||
let tool = tool_id(t).ok_or_else(|| bad(CMD, format!("unknown tool `{t}`")))?;
|
||
let mut seed = None;
|
||
let g = match tool {
|
||
"descriptive" => descriptive(s, p)?,
|
||
"histogram" => histogram(s, p)?,
|
||
"correlation" => matrix(s, p, correlation, |_| Value::Number(1.0))?,
|
||
"covariance" => matrix(s, p, covariance_p, |x| covariance_p(x, x))?,
|
||
"regression" => regression(s, p)?,
|
||
"movingAverage" => moving_average(s, p)?,
|
||
"exponentialSmoothing" => exponential_smoothing(s, p)?,
|
||
"rankPercentile" => rank_percentile(s, p)?,
|
||
"tTestPaired" | "tTestEqual" | "tTestUnequal" => t_test(s, p, tool)?,
|
||
"fTest" => f_test(s, p)?,
|
||
"anovaSingle" => anova_single(s, p)?,
|
||
"anovaTwoFactorWithout" => anova_two_factor_without(s, p)?,
|
||
"anovaTwoFactorWith" => anova_two_factor_with(s, p)?,
|
||
"zTest" => z_test(s, p)?,
|
||
"fourierAnalysis" => fourier_analysis(s, p)?,
|
||
"randomNumberGeneration" => {
|
||
// A patterned sequence draws nothing: its seed is neither read nor reported.
|
||
if !patterned(p) {
|
||
seed = Some(seed_of(p)?);
|
||
}
|
||
random_numbers(s, p, seed.unwrap_or(0))?
|
||
}
|
||
"sampling" => {
|
||
if random_sampling(p) {
|
||
seed = Some(seed_of(p)?);
|
||
}
|
||
sampling(s, p, seed)?
|
||
}
|
||
_ => return Err(bad(CMD, format!("unknown tool `{t}`"))),
|
||
};
|
||
let mut out = write_report(s, p, tool, g)?;
|
||
if let Some(k) = seed {
|
||
out["seed"] = json!(k);
|
||
}
|
||
Ok(out)
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use gridcraft_core::{CellRef, Value};
|
||
use serde_json::{Value as Json, json};
|
||
|
||
use super::{MAX_OUTCOMES, call, f64_of, pmf_table};
|
||
use crate::Session;
|
||
|
||
fn s() -> Session {
|
||
let mut s = Session::new();
|
||
s.new_workbook();
|
||
s
|
||
}
|
||
fn v(s: &Session, a: &str) -> Value {
|
||
let d = s.doc().unwrap();
|
||
let (sh, body) = crate::cmd::split_sheet(a);
|
||
let si = sh.map(|n| d.wb.sheet_index(&n).unwrap()).unwrap_or(d.wb.active_sheet());
|
||
d.wb.sheet(si).unwrap().value(CellRef::parse(body).unwrap())
|
||
}
|
||
fn n(s: &Session, a: &str) -> f64 {
|
||
v(s, a).as_f64().unwrap_or(f64::NAN)
|
||
}
|
||
fn close(a: f64, b: f64, tol: f64) -> bool {
|
||
(a - b).abs() <= tol
|
||
}
|
||
/// Evaluates a formula in a scratch cell.
|
||
fn eval(s: &mut Session, f: &str) -> f64 {
|
||
s.execute("cell.set", json!({"cell": "Z1000", "input": f})).unwrap();
|
||
let x = n(s, "Z1000");
|
||
s.execute("cell.set", json!({"cell": "Z1000", "input": ""})).unwrap();
|
||
x
|
||
}
|
||
|
||
#[test]
|
||
fn forecast_standard_errors_match_the_full_window_rescan() {
|
||
// #1043: the running sum must give what re-adding every window gave, gaps included.
|
||
let x: [f64; 11] = [2.0, 4.0, 5.0, 4.0, 5.0, 7.0, 1.0, 3.0, 8.0, 6.0, 2.0];
|
||
let f = [None, Some(3.0), Some(4.5), None, Some(4.0), Some(6.5), Some(2.0), Some(2.5), Some(7.0), Some(5.0), Some(4.0)];
|
||
let rescan = |i: usize, w: usize| -> Value {
|
||
let v = (i + 1).checked_sub(w).and_then(|s0| {
|
||
let mut acc = 0.0;
|
||
for j in s0..=i {
|
||
acc += (x.get(j)? - (*f.get(j)?)?).powi(2);
|
||
}
|
||
Some((acc / w as f64).sqrt())
|
||
});
|
||
v.map(Value::number).unwrap_or(Value::Error(gridcraft_core::CellError::NA))
|
||
};
|
||
for w in 1..=5 {
|
||
let g = super::forecast_grid(&x, &f, w, &json!({"standardErrors": true}));
|
||
for i in 0..f.len() {
|
||
let (got, want) = (g.cells.get(&(i as u32, 1)).cloned().unwrap_or(Value::Empty), rescan(i, w));
|
||
match (got.as_f64(), want.as_f64()) {
|
||
(Some(a), Some(b)) => assert!(close(a, b, 1e-12), "window {w}, row {i}: {a} vs {b}"),
|
||
_ => assert_eq!(got, want, "window {w}, row {i}"),
|
||
}
|
||
}
|
||
}
|
||
// Row 3 has no forecast: windows through it are #N/A, the next full window after it is not.
|
||
let g = super::forecast_grid(&x, &f, 3, &json!({"standardErrors": true}));
|
||
assert_eq!(g.cells.get(&(5, 1)), Some(&Value::Error(gridcraft_core::CellError::NA)));
|
||
assert!(g.cells.get(&(6, 1)).and_then(Value::as_f64).is_some());
|
||
}
|
||
|
||
#[test]
|
||
fn descriptive_statistics_of_one_to_five() {
|
||
// 1..5: mean 3, s = √2.5, kurtosis −1.2, skewness 0 (textbook values).
|
||
let mut s = s();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [["Score"], [1], [2], [3], [4], [5]]})).unwrap();
|
||
let r = s
|
||
.execute("data.analysisToolPak", json!({"tool": "Descriptive Statistics", "range": "A1:A6", "labels": true, "output": "C1", "confidence": 95, "kthLargest": 2, "kthSmallest": 1}))
|
||
.unwrap();
|
||
assert_eq!(r["range"], "C1:D18");
|
||
assert_eq!(v(&s, "C1"), Value::text("Score"));
|
||
assert_eq!(v(&s, "C3"), Value::text("Mean"));
|
||
assert_eq!(n(&s, "D3"), 3.0);
|
||
assert!(close(n(&s, "D4"), (2.5f64 / 5.0).sqrt(), 1e-12));
|
||
assert_eq!(n(&s, "D5"), 3.0);
|
||
assert_eq!(v(&s, "D6"), Value::Error(gridcraft_core::CellError::NA));
|
||
assert!(close(n(&s, "D7"), 2.5f64.sqrt(), 1e-12));
|
||
assert_eq!(n(&s, "D8"), 2.5);
|
||
assert!(close(n(&s, "D9"), -1.2, 1e-12));
|
||
assert!(close(n(&s, "D10"), 0.0, 1e-12));
|
||
assert_eq!((n(&s, "D11"), n(&s, "D12"), n(&s, "D13"), n(&s, "D14"), n(&s, "D15")), (4.0, 1.0, 5.0, 15.0, 5.0));
|
||
assert_eq!(v(&s, "C16"), Value::text("Largest(2)"));
|
||
assert_eq!(n(&s, "D16"), 4.0);
|
||
assert_eq!(n(&s, "D17"), 1.0);
|
||
assert_eq!(v(&s, "C18"), Value::text("Confidence Level(95.0%)"));
|
||
let want = eval(&mut s, "=CONFIDENCE.T(0.05,STDEV.S(A2:A6),5)");
|
||
assert!(close(n(&s, "D18"), want, 1e-9), "{} vs {want}", n(&s, "D18"));
|
||
// The output range now holds data: a second run is refused unless it may overwrite.
|
||
let again = json!({"tool": "descriptive", "range": "A2:A6", "output": "C1"});
|
||
let e = s.execute("data.analysisToolPak", again.clone()).unwrap_err().to_string();
|
||
assert!(e.contains(super::OVERWRITE_ERROR), "{e}");
|
||
let mut again = again;
|
||
again["overwrite"] = json!(true);
|
||
s.execute("data.analysisToolPak", again).unwrap();
|
||
assert_eq!(v(&s, "C1"), Value::text("Column1"));
|
||
}
|
||
|
||
#[test]
|
||
fn regression_matches_textbook_and_linest() {
|
||
// x = 1..5, y = 2,4,5,4,5: b = 0.6, a = 2.2, R² = 0.6, F = 4.5 (textbook example).
|
||
let mut s = s();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [["y", "x"], [2, 1], [4, 2], [5, 3], [4, 4], [5, 5]]})).unwrap();
|
||
s.execute(
|
||
"data.analysisToolPak",
|
||
json!({"tool": "regression", "inputY": "A1:A6", "inputX": "B1:B6", "labels": true, "newSheet": "Fit", "residuals": true, "standardizedResiduals": true}),
|
||
)
|
||
.unwrap();
|
||
assert!(close(n(&s, "Fit!B5"), 0.6, 1e-12), "R²");
|
||
assert!(close(n(&s, "Fit!B6"), 1.0 - 0.4 * 4.0 / 3.0, 1e-12), "adjusted R²");
|
||
assert!(close(n(&s, "Fit!B7"), 0.8f64.sqrt(), 1e-12), "standard error");
|
||
assert!(close(n(&s, "Fit!E12"), 4.5, 1e-12), "F");
|
||
assert!(close(n(&s, "Fit!B17"), 2.2, 1e-12), "intercept");
|
||
assert!(close(n(&s, "Fit!B18"), 0.6, 1e-12), "slope");
|
||
assert_eq!(v(&s, "Fit!A18"), Value::text("x"));
|
||
assert!(close(n(&s, "Fit!C18"), (0.8f64 / 10.0).sqrt(), 1e-12), "SE of slope");
|
||
// Simple regression: Significance F equals the slope's p-value.
|
||
assert!(close(n(&s, "Fit!F12"), n(&s, "Fit!E18"), 1e-9));
|
||
assert!(n(&s, "Fit!E18") > 0.12 && n(&s, "Fit!E18") < 0.13);
|
||
assert_eq!(v(&s, "Fit!A22"), Value::text("RESIDUAL OUTPUT"));
|
||
assert!(close(n(&s, "Fit!B25"), 2.8, 1e-12), "first prediction");
|
||
assert!(close(n(&s, "Fit!C25"), -0.8, 1e-12), "first residual");
|
||
// Multiple regression agrees with LINEST.
|
||
s.execute(
|
||
"range.setValues",
|
||
json!({"range": "Sheet1!D1", "values": [[1, 2, 3], [3, 1, 5], [4, 5, 2], [7, 3, 4], [6, 6, 1], [9, 4, 6], [8, 7, 3]]}),
|
||
)
|
||
.unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "regression", "inputY": "Sheet1!D1:D7", "inputX": "Sheet1!E1:F7", "newSheet": "Multi"}))
|
||
.unwrap();
|
||
let slope2 = eval(&mut s, "=INDEX(LINEST(Sheet1!D1:D7,Sheet1!E1:F7,TRUE,TRUE),1,1)");
|
||
let slope1 = eval(&mut s, "=INDEX(LINEST(Sheet1!D1:D7,Sheet1!E1:F7,TRUE,TRUE),1,2)");
|
||
let se1 = eval(&mut s, "=INDEX(LINEST(Sheet1!D1:D7,Sheet1!E1:F7,TRUE,TRUE),2,2)");
|
||
let icpt = eval(&mut s, "=INDEX(LINEST(Sheet1!D1:D7,Sheet1!E1:F7,TRUE,TRUE),1,3)");
|
||
assert!(close(n(&s, "Multi!B17"), icpt, 1e-9));
|
||
assert!(close(n(&s, "Multi!B18"), slope1, 1e-9));
|
||
assert!(close(n(&s, "Multi!B19"), slope2, 1e-9));
|
||
assert!(close(n(&s, "Multi!C18"), se1, 1e-9));
|
||
assert_eq!(v(&s, "Multi!A19"), Value::text("X Variable 2"));
|
||
}
|
||
|
||
#[test]
|
||
fn tests_of_means_and_variances_match_worksheet_functions() {
|
||
let mut s = s();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [[12, 10], [15, 14], [11, 9], [18, 15], [14, 13], [16, 12]]})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "tTestEqual", "range1": "A1:A6", "range2": "B1:B6", "output": "D1"})).unwrap();
|
||
assert_eq!(v(&s, "D1"), Value::text("t-Test: Two-Sample Assuming Equal Variances"));
|
||
assert_eq!(n(&s, "E9"), 10.0, "df");
|
||
let p2 = eval(&mut s, "=T.TEST(A1:A6,B1:B6,2,2)");
|
||
assert!(close(n(&s, "E13"), p2, 1e-9), "two-tail p");
|
||
s.execute("data.analysisToolPak", json!({"tool": "tTestPaired", "range1": "A1:A6", "range2": "B1:B6", "output": "H1"})).unwrap();
|
||
let pp = eval(&mut s, "=T.TEST(A1:A6,B1:B6,2,1)");
|
||
assert_eq!(v(&s, "H13"), Value::text("P(T<=t) two-tail"));
|
||
assert!(close(n(&s, "I13"), pp, 1e-9), "paired two-tail p");
|
||
s.execute("data.analysisToolPak", json!({"tool": "fTest", "range1": "A1:A6", "range2": "B1:B6", "output": "L1"})).unwrap();
|
||
let fp = eval(&mut s, "=F.TEST(A1:A6,B1:B6)");
|
||
assert!(close(n(&s, "M9") * 2.0, fp, 1e-9), "one-tail p is half of F.TEST");
|
||
// Anova on {1,2,3}, {4,5,6}, {7,8,9}: SSB 54, SSW 6, F 27.
|
||
s.execute("range.setValues", json!({"range": "A20", "values": [[1, 4, 7], [2, 5, 8], [3, 6, 9]]})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "anovaSingle", "range": "A20:C22", "output": "E20"})).unwrap();
|
||
assert_eq!(v(&s, "E31"), Value::text("Between Groups"));
|
||
assert_eq!((n(&s, "F31"), n(&s, "F32"), n(&s, "I31")), (54.0, 6.0, 27.0));
|
||
let p = eval(&mut s, "=F.DIST.RT(27,2,6)");
|
||
assert!(close(n(&s, "J31"), p, 1e-12));
|
||
}
|
||
|
||
#[test]
|
||
fn histogram_matrices_smoothing_and_ranks() {
|
||
let mut s = s();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [[1], [2], [2], [3], [5], [7]]})).unwrap();
|
||
s.execute("range.setValues", json!({"range": "B1", "values": [[2], [4], [6]]})).unwrap();
|
||
s.execute(
|
||
"data.analysisToolPak",
|
||
json!({"tool": "histogram", "range": "A1:A6", "bins": "B1:B3", "output": "D1", "cumulative": true, "pareto": true}),
|
||
)
|
||
.unwrap();
|
||
assert_eq!((n(&s, "E2"), n(&s, "E3"), n(&s, "E4"), n(&s, "E5")), (3.0, 1.0, 1.0, 1.0));
|
||
assert_eq!(v(&s, "D5"), Value::text("More"));
|
||
assert_eq!(n(&s, "F5"), 1.0);
|
||
assert_eq!((n(&s, "G2"), n(&s, "H2")), (2.0, 3.0), "Pareto puts the largest bin first");
|
||
// Correlation and covariance against CORREL and COVARIANCE.P.
|
||
s.execute("range.setValues", json!({"range": "J1", "values": [[1, 2], [2, 4], [3, 5], [4, 4], [5, 5]]})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "correlation", "range": "J1:K5", "output": "M1"})).unwrap();
|
||
assert_eq!(n(&s, "N2"), 1.0);
|
||
assert!(close(n(&s, "N3"), eval(&mut s, "=CORREL(J1:J5,K1:K5)"), 1e-12));
|
||
assert!(v(&s, "O2").is_empty(), "upper triangle is blank");
|
||
s.execute("data.analysisToolPak", json!({"tool": "covariance", "range": "J1:K5", "output": "M10"})).unwrap();
|
||
assert!(close(n(&s, "N12"), eval(&mut s, "=COVARIANCE.P(J1:J5,K1:K5)"), 1e-12));
|
||
assert!(close(n(&s, "N11"), eval(&mut s, "=VAR.P(J1:J5)"), 1e-12));
|
||
// Moving average of 3 and exponential smoothing with damping 0.3.
|
||
s.execute("data.analysisToolPak", json!({"tool": "movingAverage", "range": "K1:K5", "output": "Q1", "standardErrors": true})).unwrap();
|
||
assert_eq!(v(&s, "Q2"), Value::Error(gridcraft_core::CellError::NA));
|
||
assert!(close(n(&s, "Q3"), 11.0 / 3.0, 1e-12));
|
||
assert!(close(n(&s, "Q5"), 14.0 / 3.0, 1e-12));
|
||
assert_eq!(v(&s, "R4"), Value::Error(gridcraft_core::CellError::NA));
|
||
assert!(n(&s, "R5").is_finite());
|
||
s.execute("data.analysisToolPak", json!({"tool": "exponentialSmoothing", "range": "K1:K5", "output": "S1"})).unwrap();
|
||
assert_eq!(n(&s, "S2"), 2.0);
|
||
assert!(close(n(&s, "S3"), 0.7 * 4.0 + 0.3 * 2.0, 1e-12));
|
||
// Rank and Percentile of 2,4,5,4,5: the two 5s share rank 1 (3 of 4 others below them).
|
||
s.execute("data.analysisToolPak", json!({"tool": "rankPercentile", "range": "K1:K5", "output": "U1"})).unwrap();
|
||
assert_eq!((n(&s, "U2"), n(&s, "V2"), n(&s, "W2"), n(&s, "X2")), (3.0, 5.0, 1.0, 0.75));
|
||
assert_eq!((n(&s, "U3"), n(&s, "W3")), (5.0, 1.0));
|
||
assert_eq!((n(&s, "V6"), n(&s, "X6")), (2.0, 0.0));
|
||
assert_eq!(n(&s, "X4"), 0.25, "4 has one value below it among five");
|
||
}
|
||
|
||
#[test]
|
||
fn two_factor_anova_without_replication_matches_a_textbook_example() {
|
||
// Hardness testing, a randomized complete block design (Montgomery, Design and Analysis
|
||
// of Experiments, Example 4.1): four tips (rows) by four coupons (columns). Published:
|
||
// SS tips 0.385, SS coupons 0.825, SS error 0.080, SS total 1.290, F tips 14.44.
|
||
let mut s = s();
|
||
s.execute(
|
||
"range.setValues",
|
||
json!({"range": "A1", "values": [
|
||
["", "Coupon 1", "Coupon 2", "Coupon 3", "Coupon 4"],
|
||
["Tip 1", 9.3, 9.4, 9.6, 10.0],
|
||
["Tip 2", 9.4, 9.3, 9.8, 9.9],
|
||
["Tip 3", 9.2, 9.4, 9.5, 9.7],
|
||
["Tip 4", 9.7, 9.6, 10.0, 10.2],
|
||
]}),
|
||
)
|
||
.unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "Anova: Two-Factor Without Replication", "range": "A1:E5", "labels": true, "output": "H1"}))
|
||
.unwrap();
|
||
assert_eq!(v(&s, "H1"), Value::text("Anova: Two-Factor Without Replication"));
|
||
assert_eq!((v(&s, "H4"), n(&s, "I4")), (Value::text("Tip 1"), 4.0));
|
||
assert!(close(n(&s, "J4"), 38.3, 1e-12), "sum of tip 1");
|
||
assert_eq!(v(&s, "H9"), Value::text("Coupon 1"));
|
||
assert!(close(n(&s, "K12"), 9.95, 1e-12), "average of coupon 4");
|
||
assert_eq!(v(&s, "H17"), Value::text("Rows"));
|
||
assert!(close(n(&s, "I17"), 0.385, 1e-12));
|
||
assert!(close(n(&s, "I18"), 0.825, 1e-12));
|
||
assert!(close(n(&s, "I19"), 0.080, 1e-12));
|
||
assert!(close(n(&s, "I21"), 1.290, 1e-12));
|
||
assert_eq!((n(&s, "J17"), n(&s, "J19"), n(&s, "J21")), (3.0, 9.0, 15.0));
|
||
assert!(close(n(&s, "L17"), 14.4375, 1e-9), "F for tips");
|
||
assert!(close(n(&s, "M17"), eval(&mut s, "=F.DIST.RT(14.4375,3,9)"), 1e-9));
|
||
assert!(close(n(&s, "N17"), eval(&mut s, "=F.INV.RT(0.05,3,9)"), 1e-9));
|
||
assert!(v(&s, "L19").is_empty(), "the error line has no F");
|
||
}
|
||
|
||
#[test]
|
||
fn two_factor_anova_with_replication_matches_a_textbook_example() {
|
||
let mut s = s();
|
||
// Battery life (Montgomery, Example 5.1): three materials (samples of four rows) by three
|
||
// temperatures (columns). Published: SS material 10683.72, SS temperature 39118.72,
|
||
// SS interaction 9613.78, SS error 18230.75 on 27 df, SS total 77646.97.
|
||
s.execute(
|
||
"range.setValues",
|
||
json!({"range": "A20", "values": [
|
||
["Material", "15 F", "70 F", "125 F"],
|
||
["Type 1", 130, 34, 20], ["", 155, 40, 70], ["", 74, 80, 82], ["", 180, 75, 58],
|
||
["Type 2", 150, 136, 25], ["", 188, 122, 70], ["", 159, 106, 58], ["", 126, 115, 45],
|
||
["Type 3", 138, 174, 96], ["", 110, 120, 104], ["", 168, 150, 82], ["", 160, 139, 60],
|
||
]}),
|
||
)
|
||
.unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "anovaTwoFactorWith", "range": "A20:D32", "rowsPerSample": 4, "newSheet": "Battery"}))
|
||
.unwrap();
|
||
assert_eq!(v(&s, "Battery!A3"), Value::text("SUMMARY"));
|
||
assert_eq!((v(&s, "Battery!B3"), v(&s, "Battery!E3")), (Value::text("15 F"), Value::text("Total")));
|
||
assert_eq!(v(&s, "Battery!A4"), Value::text("Type 1"));
|
||
assert_eq!((n(&s, "Battery!B5"), n(&s, "Battery!B6"), n(&s, "Battery!E5")), (4.0, 539.0, 12.0));
|
||
assert_eq!(v(&s, "Battery!A22"), Value::text("Total"));
|
||
assert_eq!(n(&s, "Battery!B24"), 1738.0, "sum of the 15 F column");
|
||
assert_eq!(v(&s, "Battery!A29"), Value::text("ANOVA"));
|
||
assert_eq!(
|
||
(v(&s, "Battery!A31"), v(&s, "Battery!A34"), v(&s, "Battery!A36")),
|
||
(Value::text("Sample"), Value::text("Within"), Value::text("Total"))
|
||
);
|
||
assert!(close(n(&s, "Battery!B31"), 10683.72, 0.01));
|
||
assert!(close(n(&s, "Battery!B32"), 39118.72, 0.01));
|
||
assert!(close(n(&s, "Battery!B33"), 9613.78, 0.01));
|
||
assert!(close(n(&s, "Battery!B34"), 18230.75, 0.01));
|
||
assert!(close(n(&s, "Battery!B36"), 77646.97, 0.01));
|
||
assert_eq!((n(&s, "Battery!C33"), n(&s, "Battery!C34"), n(&s, "Battery!C36")), (4.0, 27.0, 35.0));
|
||
assert!(close(n(&s, "Battery!E31"), 7.91, 0.005) && close(n(&s, "Battery!E32"), 28.97, 0.005) && close(n(&s, "Battery!E33"), 3.56, 0.005));
|
||
assert!(close(n(&s, "Battery!F33"), 0.0186, 0.00005), "interaction P-value");
|
||
// Rows per sample must divide the data rows.
|
||
let bad = json!({"tool": "anovaTwoFactorWith", "range": "Sheet1!A20:D32", "rowsPerSample": 5, "output": "F60"});
|
||
assert!(s.execute("data.analysisToolPak", bad).is_err());
|
||
}
|
||
|
||
#[test]
|
||
fn z_test_with_known_variances() {
|
||
// x̄ = 3 (n = 5, σ² = 2), ȳ = 4 (n = 3, σ² = 6): z = −1/√2.4. Critical values from the
|
||
// standard normal table: 1.644854 (one tail) and 1.959964 (two tails) at α = 0.05.
|
||
let mut s = s();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [["x", "y"], [1, 2], [2, 4], [3, 6], [4, null], [5, null]]})).unwrap();
|
||
s.execute(
|
||
"data.analysisToolPak",
|
||
json!({"tool": "zTest", "range1": "A1:A6", "range2": "B1:B6", "labels": true, "variance1": 2, "variance2": 6, "output": "D1"}),
|
||
)
|
||
.unwrap();
|
||
assert_eq!(v(&s, "D1"), Value::text("z-Test: Two Sample for Means"));
|
||
assert_eq!((v(&s, "E3"), v(&s, "F3")), (Value::text("x"), Value::text("y")));
|
||
assert_eq!((n(&s, "E4"), n(&s, "F4"), n(&s, "E5"), n(&s, "F6")), (3.0, 4.0, 2.0, 3.0));
|
||
let z = -1.0 / 2.4f64.sqrt();
|
||
assert!(close(n(&s, "E8"), z, 1e-12));
|
||
let tail = 1.0 - eval(&mut s, &format!("=NORM.S.DIST({},TRUE)", -z));
|
||
assert!(close(n(&s, "E9"), tail, 1e-12) && close(n(&s, "E9"), 0.2593, 1e-4));
|
||
assert!(close(n(&s, "E10"), 1.644_853_6, 1e-6));
|
||
assert!(close(n(&s, "E11"), 2.0 * tail, 1e-12));
|
||
assert!(close(n(&s, "E12"), 1.959_964_0, 1e-6));
|
||
// Unknown or non-positive variances are refused.
|
||
assert!(s.execute("data.analysisToolPak", json!({"tool": "zTest", "range1": "A2:A6", "range2": "B2:B4", "output": "H1"})).is_err());
|
||
let zero = json!({"tool": "zTest", "range1": "A2:A6", "range2": "B2:B4", "variance1": 0, "variance2": 1, "output": "H1"});
|
||
assert!(s.execute("data.analysisToolPak", zero).is_err());
|
||
}
|
||
|
||
#[test]
|
||
fn fourier_analysis_matches_known_transforms() {
|
||
let mut s = s();
|
||
// DFT of 1, 2, 3, 4: 10, −2+2i, −2, −2−2i.
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [["Signal"], [1], [2], [3], [4]]})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "Fourier Analysis", "range": "A1:A5", "labels": true, "output": "C1"})).unwrap();
|
||
let col = |s: &Session, c: &str, rows: std::ops::RangeInclusive<u32>| rows.map(|r| v(s, &format!("{c}{r}"))).collect::<Vec<_>>();
|
||
assert_eq!(col(&s, "C", 1..=5), ["Signal", "10", "-2+2i", "-2", "-2-2i"].map(Value::text).to_vec());
|
||
// The inverse transform of that output returns the signal.
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "C2:C5", "inverse": true, "output": "E1"})).unwrap();
|
||
assert_eq!(col(&s, "E", 1..=4), ["1", "2", "3", "4"].map(Value::text).to_vec());
|
||
// cos(2πn/8) over eight points: N/2 = 4 at the first and last frequencies, 0 elsewhere.
|
||
let wave: Vec<Json> = (0..8).map(|k| json!([(std::f64::consts::PI * k as f64 / 4.0).cos()])).collect();
|
||
s.execute("range.setValues", json!({"range": "G1", "values": wave})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "G1:G8", "output": "H1"})).unwrap();
|
||
assert_eq!(col(&s, "H", 1..=8), ["0", "4", "0", "0", "0", "0", "0", "4"].map(Value::text).to_vec());
|
||
// Complex text input; a single value is its own transform.
|
||
s.execute("range.setValues", json!({"range": "J1", "values": [["3+4i"], ["1-i"]]})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "J1:J2", "output": "K1"})).unwrap();
|
||
assert_eq!(col(&s, "K", 1..=2), ["4+3i", "2+5i"].map(Value::text).to_vec());
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "J1", "output": "L1"})).unwrap();
|
||
assert_eq!(v(&s, "L1"), Value::text("3+4i"));
|
||
// Three values are not a power of two; a word is not a number.
|
||
assert!(s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "A2:A4", "output": "N1"})).is_err());
|
||
assert!(s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "A1:A4", "output": "N1"})).is_err());
|
||
assert!(s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "A1:A8192", "output": "N1"})).is_err());
|
||
}
|
||
|
||
#[test]
|
||
fn fourier_analysis_keeps_a_small_imaginary_part_beside_a_large_real_one() {
|
||
let mut s = s();
|
||
// With one value the transform is the identity and performs no arithmetic, so no part of
|
||
// the output is rounding noise: 1e14 + 1i comes back as 1e14 + 1i, forward and inverse.
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [["100000000000000+1i"]]})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "A1", "output": "C1"})).unwrap();
|
||
assert_eq!(eval(&mut s, "=IMREAL(C1)"), 1e14);
|
||
assert_eq!(eval(&mut s, "=IMAGINARY(C1)"), 1.0);
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "C1", "inverse": true, "output": "E1"})).unwrap();
|
||
assert_eq!(eval(&mut s, "=IMREAL(E1)"), 1e14);
|
||
assert_eq!(eval(&mut s, "=IMAGINARY(E1)"), 1.0);
|
||
}
|
||
|
||
#[test]
|
||
fn fourier_analysis_round_trip_recovers_small_components_beside_a_large_dc() {
|
||
let mut s = s();
|
||
// x_n = 1e6 + 1e-7·cos(2πn/8) for n = 0 … 7. The transform of a constant c over N points is
|
||
// N·c at k = 0 and 0 elsewhere; the transform of a·cos(2πn/N) is a·N/2 at k = 1 and
|
||
// k = N − 1. So X_0 = 8·1e6 = 8e6, X_1 = X_7 = 1e-7·8/2 = 4e-7 and X_2 … X_6 = 0. The inputs
|
||
// are stored to within ulp(1e6)/2 ≈ 6e-11 and the first butterfly stage subtracts nearby
|
||
// numbers exactly, so the harmonics carry at most about 5e-10 of noise. The cleanup bound
|
||
// 8·ε·log2(8)·Σ|x_n| ≈ 8 · 2.2e-16 · 3 · 8e6 ≈ 4.3e-8 zeroes noise of that size and keeps
|
||
// 4e-7; a cutoff relative to the largest output (8e6 · 1e-13 = 8e-7) would zero the
|
||
// harmonics.
|
||
let x = |k: usize| 1e6 + 1e-7 * (std::f64::consts::PI * k as f64 / 4.0).cos();
|
||
let wave: Vec<Json> = (0..8).map(|k| json!([x(k)])).collect();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": wave})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "A1:A8", "output": "C1"})).unwrap();
|
||
assert!(close(eval(&mut s, "=IMREAL(C1)"), 8e6, 1e-6));
|
||
assert_eq!(eval(&mut s, "=IMAGINARY(C1)"), 0.0);
|
||
for cell in ["C2", "C8"] {
|
||
assert!(close(eval(&mut s, &format!("=IMREAL({cell})")), 4e-7, 5e-9), "{cell}");
|
||
assert_eq!(eval(&mut s, &format!("=IMAGINARY({cell})")), 0.0, "{cell}");
|
||
}
|
||
for r in 3..=7 {
|
||
assert_eq!(v(&s, &format!("C{r}")), Value::text("0"), "C{r}");
|
||
}
|
||
// The inverse of that output returns each x_n to within 2e-9 (2e-15 relative); without the
|
||
// harmonics every value would come back as 1e6, off by up to 1e-7.
|
||
s.execute("data.analysisToolPak", json!({"tool": "fourierAnalysis", "range": "C1:C8", "inverse": true, "output": "E1"})).unwrap();
|
||
for k in 0..8 {
|
||
let cell = format!("E{}", k + 1);
|
||
assert!(close(eval(&mut s, &format!("=IMREAL({cell})")), x(k), 2e-9), "{cell}");
|
||
assert_eq!(eval(&mut s, &format!("=IMAGINARY({cell})")), 0.0, "{cell}");
|
||
}
|
||
}
|
||
|
||
/// Column `col` of sheet `sheet`, rows 1 to `rows`, as numbers.
|
||
fn numbers_of(s: &Session, sheet: &str, col: &str, rows: u32) -> Vec<f64> {
|
||
(1..=rows).map(|r| n(s, &format!("{sheet}!{col}{r}"))).collect()
|
||
}
|
||
fn mean_var(x: &[f64]) -> (f64, f64) {
|
||
let m = x.iter().sum::<f64>() / x.len() as f64;
|
||
(m, x.iter().map(|v| (v - m).powi(2)).sum::<f64>() / (x.len() - 1) as f64)
|
||
}
|
||
|
||
#[test]
|
||
fn random_numbers_are_seeded_and_follow_their_distributions() {
|
||
let mut s = s();
|
||
let rng = |s: &mut Session, sheet: &str, p: Json| {
|
||
let mut p = p;
|
||
p["tool"] = json!("randomNumberGeneration");
|
||
p["newSheet"] = json!(sheet);
|
||
s.execute("data.analysisToolPak", p).unwrap()
|
||
};
|
||
// The same seed repeats the numbers; the result reports the seed.
|
||
let r = rng(&mut s, "U1", json!({"variables": 2, "count": 2000, "distribution": "uniform", "from": 2, "to": 4, "seed": 42}));
|
||
assert_eq!((r["seed"].as_u64(), r["range"].as_str()), (Some(42), Some("A1:B2000")));
|
||
rng(&mut s, "U2", json!({"variables": 2, "count": 2000, "distribution": "uniform", "from": 2, "to": 4, "seed": 42}));
|
||
let (a, b) = (numbers_of(&s, "U1", "B", 2000), numbers_of(&s, "U2", "B", 2000));
|
||
assert_eq!(a, b);
|
||
assert!(a.iter().all(|x| (2.0..4.0).contains(x)));
|
||
// Uniform on [2, 4]: mean 3, variance 1/3.
|
||
let (m, var) = mean_var(&a);
|
||
assert!((m - 3.0).abs() < 0.05 && (var - 1.0 / 3.0).abs() < 0.03, "{m} {var}");
|
||
let unseeded = rng(&mut s, "U3", json!({"count": 1}));
|
||
assert!(unseeded["seed"].as_u64().is_some(), "a seed is chosen and reported");
|
||
// Normal(10, 2): mean 10, variance 4.
|
||
rng(&mut s, "N", json!({"count": 20000, "distribution": "normal", "mean": 10, "standardDeviation": 2, "seed": 7}));
|
||
let (m, var) = mean_var(&numbers_of(&s, "N", "A", 20000));
|
||
assert!((m - 10.0).abs() < 0.05 && (var - 4.0).abs() < 0.15, "{m} {var}");
|
||
// Binomial(20, 0.3): mean np = 6, variance np(1 − p) = 4.2.
|
||
rng(&mut s, "B", json!({"count": 20000, "distribution": "binomial", "p": 0.3, "trials": 20, "seed": 8}));
|
||
let x = numbers_of(&s, "B", "A", 20000);
|
||
assert!(x.iter().all(|v| v.fract() == 0.0 && (0.0..=20.0).contains(v)));
|
||
let (m, var) = mean_var(&x);
|
||
assert!((m - 6.0).abs() < 0.08 && (var - 4.2).abs() < 0.2, "{m} {var}");
|
||
// Poisson(3.5): mean and variance 3.5.
|
||
rng(&mut s, "P", json!({"count": 20000, "distribution": "poisson", "lambda": 3.5, "seed": 9}));
|
||
let (m, var) = mean_var(&numbers_of(&s, "P", "A", 20000));
|
||
assert!((m - 3.5).abs() < 0.08 && (var - 3.5).abs() < 0.2, "{m} {var}");
|
||
// Bernoulli(0.25) and a discrete distribution on {1, 2} with P(2) = 0.75.
|
||
rng(&mut s, "Ber", json!({"count": 4000, "distribution": "bernoulli", "p": 0.25, "seed": 10}));
|
||
let ones = numbers_of(&s, "Ber", "A", 4000).iter().filter(|v| **v == 1.0).count();
|
||
assert!((900..1100).contains(&ones), "{ones}");
|
||
s.execute("range.setValues", json!({"range": "Sheet1!A1", "values": [[1, 0.25], [2, 0.75]]})).unwrap();
|
||
rng(&mut s, "D", json!({"count": 4000, "distribution": "discrete", "valueRange": "Sheet1!A1:B2", "seed": 11}));
|
||
let x = numbers_of(&s, "D", "A", 4000);
|
||
assert!(x.iter().all(|v| *v == 1.0 || *v == 2.0));
|
||
let twos = x.iter().filter(|v| **v == 2.0).count();
|
||
assert!((2900..3100).contains(&twos), "{twos}");
|
||
// Patterned: 1 to 3, each number twice, the sequence twice.
|
||
rng(&mut s, "Pat", json!({"distribution": "patterned", "from": 1, "to": 3, "step": 1, "repeatEach": 2, "repeatSequence": 2}));
|
||
assert_eq!(numbers_of(&s, "Pat", "A", 12), [1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 1.0, 1.0, 2.0, 2.0, 3.0, 3.0]);
|
||
assert!(v(&s, "Pat!A13").is_empty());
|
||
// Output is bounded; parameters are checked.
|
||
for p in [
|
||
json!({"count": 2_000_000}),
|
||
json!({"count": 1000, "variables": 2000}),
|
||
json!({"count": 0}),
|
||
json!({"count": 1.5}),
|
||
json!({"distribution": "uniform"}),
|
||
json!({"count": 1, "distribution": "bernoulli", "p": 2}),
|
||
json!({"count": 1, "distribution": "binomial", "trials": 1e9}),
|
||
json!({"count": 1, "distribution": "poisson", "lambda": -1}),
|
||
json!({"count": 1, "distribution": "normal", "standardDeviation": -1}),
|
||
json!({"count": 1, "distribution": "discrete", "valueRange": "Sheet1!A1:A2"}),
|
||
json!({"count": 1, "distribution": "patterned", "from": 1, "to": 1e9, "step": 1}),
|
||
json!({"count": 1, "distribution": "patterned", "from": 3, "to": 1}),
|
||
json!({"count": 1, "distribution": "weird"}),
|
||
json!({"count": 1, "seed": -1}),
|
||
json!({"count": 1, "seed": "abc"}),
|
||
] {
|
||
let mut q = p.clone();
|
||
q["tool"] = json!("randomNumberGeneration");
|
||
q["output"] = json!("Sheet1!H1");
|
||
assert!(s.execute("data.analysisToolPak", q).is_err(), "{p}");
|
||
}
|
||
s.execute("range.setValues", json!({"range": "Sheet1!A1", "values": [[1, 0.5], [2, 0.6]]})).unwrap();
|
||
let p =
|
||
json!({"tool": "randomNumberGeneration", "count": 1, "distribution": "discrete", "valueRange": "Sheet1!A1:B2", "output": "Sheet1!H1"});
|
||
assert!(s.execute("data.analysisToolPak", p).is_err(), "probabilities must add up to 1");
|
||
}
|
||
|
||
#[test]
|
||
fn binomial_and_poisson_tables_agree_with_the_worksheet_functions() {
|
||
// The recurrence outward from the mode reproduces BINOM.DIST and POISSON.DIST at every
|
||
// tabulated outcome, the table carries the whole distribution, and it is only as wide as
|
||
// the outcomes that matter. The worksheet functions themselves are accurate to about
|
||
// n · ε for n trials (their logarithms are of that size), hence the moderate n and λ.
|
||
let binom = |k: f64, n: f64, p: f64| f64_of(&call("BINOM.DIST", &[Value::Number(k), Value::Number(n), Value::Number(p), Value::Bool(false)]));
|
||
for (n, p) in [(20.0f64, 0.3f64), (100.0, 0.5), (50.0, 0.9), (30.0, 0.01), (7.0, 0.999), (1.0, 0.5)] {
|
||
let mode = ((n + 1.0) * p).floor().min(n);
|
||
let odds = p / (1.0 - p);
|
||
let (first, pmf) = pmf_table(n, mode, binom(mode, n, p), |k| (n - k) / (k + 1.0) * odds).unwrap();
|
||
assert!(first >= 0.0 && first + pmf.len() as f64 - 1.0 <= n, "n={n} p={p}: {first} + {}", pmf.len());
|
||
for (i, got) in pmf.iter().enumerate() {
|
||
let k = first + i as f64;
|
||
let want = binom(k, n, p);
|
||
assert!(((got - want) / want).abs() < 1e-12, "n={n} p={p} k={k}: {got} vs {want}");
|
||
}
|
||
assert!((pmf.iter().sum::<f64>() - 1.0).abs() < 1e-12, "n={n} p={p}");
|
||
}
|
||
let poisson = |k: f64, lambda: f64| f64_of(&call("POISSON.DIST", &[Value::Number(k), Value::Number(lambda), Value::Bool(false)]));
|
||
for lambda in [0.2f64, 3.5, 50.0] {
|
||
let mode = lambda.floor();
|
||
let (first, pmf) = pmf_table(MAX_OUTCOMES, mode, poisson(mode, lambda), |k| lambda / (k + 1.0)).unwrap();
|
||
assert_eq!(first, 0.0);
|
||
for (i, got) in pmf.iter().enumerate() {
|
||
let want = poisson(i as f64, lambda);
|
||
assert!(((got - want) / want).abs() < 1e-12, "lambda={lambda} k={i}: {got} vs {want}");
|
||
}
|
||
assert!((pmf.iter().sum::<f64>() - 1.0).abs() < 1e-12, "lambda={lambda}");
|
||
}
|
||
// Binomial(1e6, 0.5): standard deviation 500, so the terms vanish within about 38
|
||
// standard deviations below the mode and the sum reaches 1 within about 8 above it:
|
||
// roughly 23,000 entries starting near 481,000, not 1,000,001 starting at 0.
|
||
let (n, p) = (1_000_000.0, 0.5);
|
||
let (first, pmf) = pmf_table(n, 500_000.0, binom(500_000.0, n, p), |k| (n - k) / (k + 1.0)).unwrap();
|
||
assert!((470_000.0..=490_000.0).contains(&first) && pmf.len() < 40_000, "{first} + {}", pmf.len());
|
||
// The anchor BINOM.DIST(500000, 1e6, 0.5) carries that function's error of about n · ε
|
||
// (its logarithms are near 7e5), which scales the whole table; draws divide by the total.
|
||
assert!((pmf.iter().sum::<f64>() - 1.0).abs() < 1e-8);
|
||
// A bad anchor or mode is refused rather than tabulated.
|
||
assert!(pmf_table(10.0, 5.0, f64::NAN, |_| 1.0).is_err());
|
||
assert!(pmf_table(10.0, 5.0, 0.0, |_| 1.0).is_err());
|
||
assert!(pmf_table(10.0, 11.0, 0.5, |_| 1.0).is_err());
|
||
}
|
||
|
||
#[test]
|
||
#[allow(clippy::disallowed_methods)] // a native-only timing test: the clock is read off wasm
|
||
fn binomial_and_poisson_draws_at_the_limits_are_tabulated_quickly() {
|
||
let mut s = s();
|
||
// The tables cover the outcomes around the mode, not every one of the 1,000,001 for a
|
||
// million trials, so the largest inputs cost milliseconds. Binomial(1e6, 0.5) has mean
|
||
// 500000 and standard deviation 500: 1000 draws put the sample mean within 500/√1000 ≈ 16
|
||
// of it, so 100 is six standard errors. Poisson(1e5) has mean 1e5 and standard deviation
|
||
// 316: 1000 draws put the mean within 10, so 60 is six standard errors. With p = 1 every
|
||
// draw is exactly 1e6, with p = 0 exactly 0.
|
||
let mut draw = |sheet: &str, dist: Json| {
|
||
let mut dist = dist;
|
||
dist["tool"] = json!("randomNumberGeneration");
|
||
dist["count"] = json!(1000);
|
||
dist["seed"] = json!(3);
|
||
dist["newSheet"] = json!(sheet);
|
||
s.execute("data.analysisToolPak", dist).unwrap();
|
||
};
|
||
let start = std::time::Instant::now();
|
||
draw("Half", json!({"distribution": "binomial", "p": 0.5, "trials": 1_000_000}));
|
||
draw("Sure", json!({"distribution": "binomial", "p": 1, "trials": 1_000_000}));
|
||
draw("Never", json!({"distribution": "binomial", "p": 0, "trials": 1_000_000}));
|
||
draw("Pois", json!({"distribution": "poisson", "lambda": 100_000}));
|
||
let elapsed = start.elapsed();
|
||
let (m, _) = mean_var(&numbers_of(&s, "Half", "A", 1000));
|
||
assert!((m - 500_000.0).abs() < 100.0, "{m}");
|
||
assert!(numbers_of(&s, "Sure", "A", 1000).iter().all(|v| *v == 1_000_000.0));
|
||
assert!(numbers_of(&s, "Never", "A", 1000).iter().all(|v| *v == 0.0));
|
||
let (m, _) = mean_var(&numbers_of(&s, "Pois", "A", 1000));
|
||
assert!((m - 100_000.0).abs() < 60.0, "{m}");
|
||
assert!(elapsed < std::time::Duration::from_secs(2), "{elapsed:?}");
|
||
}
|
||
|
||
#[test]
|
||
fn sampling_periodic_and_random() {
|
||
let mut s = s();
|
||
let values: Vec<Json> = (1..=10).map(|k| json!([k])).collect();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": values})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "Sampling", "range": "A1:A10", "period": 3, "output": "C1"})).unwrap();
|
||
assert_eq!(numbers_of(&s, "Sheet1", "C", 3), [3.0, 6.0, 9.0]);
|
||
assert!(v(&s, "C4").is_empty());
|
||
let r = s
|
||
.execute(
|
||
"data.analysisToolPak",
|
||
json!({"tool": "sampling", "range": "A1:A10", "method": "random", "samples": 20, "seed": 5, "output": "D1"}),
|
||
)
|
||
.unwrap();
|
||
assert_eq!(r["seed"], 5);
|
||
s.execute(
|
||
"data.analysisToolPak",
|
||
json!({"tool": "sampling", "range": "A1:A10", "method": "random", "samples": 20, "seed": 5, "output": "E1"}),
|
||
)
|
||
.unwrap();
|
||
let (a, b) = (numbers_of(&s, "Sheet1", "D", 20), numbers_of(&s, "Sheet1", "E", 20));
|
||
assert_eq!(a, b, "the same seed repeats the sample");
|
||
assert!(a.iter().all(|x| x.fract() == 0.0 && (1.0..=10.0).contains(x)));
|
||
for p in [
|
||
json!({"range": "A1:A10", "period": 11}),
|
||
json!({"range": "A1:A10", "period": 0}),
|
||
json!({"range": "A1:A10"}),
|
||
json!({"range": "A1:A10", "method": "random", "samples": 2_000_000}),
|
||
json!({"range": "B1:B10", "period": 1}),
|
||
] {
|
||
let mut q = p.clone();
|
||
q["tool"] = json!("sampling");
|
||
q["output"] = json!("G1");
|
||
assert!(s.execute("data.analysisToolPak", q).is_err(), "{p}");
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn options_typed_as_text_that_is_not_a_number_are_refused() {
|
||
// A form sends its options as text; one that doesn't read as a number (a decimal comma,
|
||
// a typo) is an error, not the default, while a blank option takes the default.
|
||
let mut s = s();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [[1, 2], [2, 4], [3, 6]]})).unwrap();
|
||
let rng = |dist: &str| json!({"tool": "randomNumberGeneration", "count": 3, "distribution": dist, "seed": 1, "outputTo": "newSheet"});
|
||
let z = json!({"tool": "zTest", "range1": "Sheet1!A1:A3", "range2": "Sheet1!B1:B3", "variance1": 1, "variance2": 4, "outputTo": "newSheet"});
|
||
let with = |base: &Json, key: &str, val: &str| {
|
||
let mut p = base.clone();
|
||
p[key] = json!(val);
|
||
p
|
||
};
|
||
let cases = [
|
||
(rng("uniform"), "from"),
|
||
(rng("uniform"), "to"),
|
||
(rng("normal"), "mean"),
|
||
(rng("normal"), "standardDeviation"),
|
||
(rng("bernoulli"), "p"),
|
||
(rng("binomial"), "p"),
|
||
(rng("poisson"), "lambda"),
|
||
(rng("patterned"), "step"),
|
||
(rng("uniform"), "variables"),
|
||
(z.clone(), "hypothesizedDifference"),
|
||
];
|
||
for (base, key) in cases {
|
||
let e = s.execute("data.analysisToolPak", with(&base, key, "0,3")).unwrap_err().to_string();
|
||
assert!(e.contains(key), "{key}: {e}");
|
||
s.execute("data.analysisToolPak", with(&base, key, " ")).unwrap_or_else(|e| panic!("{key} blank: {e}"));
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn an_alpha_that_is_not_a_number_is_refused() {
|
||
// A decimal comma or a typo in Alpha must not run the test at 0.05 without a word.
|
||
let mut s = s();
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [["", "c1", "c2"], ["r1", 1, 2], ["r1", 2, 4], ["r2", 3, 6], ["r2", 5, 9]]}))
|
||
.unwrap();
|
||
for base in [
|
||
json!({"tool": "zTest", "range1": "Sheet1!B2:B5", "range2": "Sheet1!C2:C5", "variance1": 1, "variance2": 4}),
|
||
json!({"tool": "anovaTwoFactorWithout", "range": "Sheet1!B2:C5"}),
|
||
json!({"tool": "anovaTwoFactorWith", "range": "Sheet1!A1:C5", "rowsPerSample": 2}),
|
||
json!({"tool": "tTestEqual", "range1": "Sheet1!B2:B5", "range2": "Sheet1!C2:C5"}),
|
||
] {
|
||
let mut p = base.clone();
|
||
p["alpha"] = json!("0,01");
|
||
let e = s.execute("data.analysisToolPak", p).unwrap_err().to_string();
|
||
assert!(e.contains("alpha"), "{base}: {e}");
|
||
let mut p = base.clone();
|
||
p["alpha"] = json!("");
|
||
s.execute("data.analysisToolPak", p).unwrap_or_else(|e| panic!("{base} blank alpha: {e}"));
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn a_patterned_sequence_ignores_the_seed() {
|
||
// Patterned generation draws nothing, so a seed left in the form (even one that isn't a
|
||
// number) neither blocks the run nor shows up in the result.
|
||
let mut s = s();
|
||
for seed in [json!("abc"), json!(7)] {
|
||
let r = s
|
||
.execute(
|
||
"data.analysisToolPak",
|
||
json!({"tool": "randomNumberGeneration", "distribution": "patterned", "from": 1, "to": 3, "seed": seed, "outputTo": "newSheet"}),
|
||
)
|
||
.unwrap_or_else(|e| panic!("seed {seed}: {e}"));
|
||
assert!(r.get("seed").is_none(), "seed {seed}: {r}");
|
||
}
|
||
// A random distribution still refuses the bad seed.
|
||
let bad = json!({"tool": "randomNumberGeneration", "distribution": "uniform", "count": 2, "seed": "abc", "outputTo": "newSheet"});
|
||
assert!(s.execute("data.analysisToolPak", bad).is_err());
|
||
}
|
||
|
||
#[test]
|
||
fn hostile_inputs_are_refused_without_panicking() {
|
||
let mut s = s();
|
||
let r = s.execute("data.analysisToolPak", json!({})).unwrap();
|
||
assert!(r["tools"].as_array().is_some_and(|t| t.len() == super::TOOLS.len()));
|
||
assert!(s.take_ui_requests().iter().any(|q| matches!(q, crate::UiRequest::Dialog(n, _) if n == "analysisToolPak")));
|
||
s.execute("range.setValues", json!({"range": "A1", "values": [[1, "x"], [2, 3]]})).unwrap();
|
||
for p in [
|
||
json!({"tool": "nope"}),
|
||
json!({"tool": "descriptive"}),
|
||
json!({"tool": "descriptive", "range": "A1:XFD1048576"}),
|
||
json!({"tool": "descriptive", "range": "A1:B2"}),
|
||
json!({"tool": "descriptive", "range": "A1:A2", "output": "XFD1048576"}),
|
||
json!({"tool": "descriptive", "range": "A1:A2", "newSheet": "bad[name]"}),
|
||
json!({"tool": "regression", "inputY": "A1:A2", "inputX": "A1:A2"}),
|
||
json!({"tool": "regression", "inputY": "A1:A2", "inputX": "A1:Q2"}),
|
||
json!({"tool": "movingAverage", "range": "A1:A2", "interval": 1e300}),
|
||
json!({"tool": "movingAverage", "range": "A1:A2", "interval": -3}),
|
||
json!({"tool": "exponentialSmoothing", "range": "A1:A2", "damping": 7}),
|
||
json!({"tool": "histogram", "range": "C1:C3"}),
|
||
json!({"tool": "tTestPaired", "range1": "A1", "range2": "A2"}),
|
||
json!({"tool": "anovaSingle", "range": "A1:A2"}),
|
||
json!({"tool": "correlation", "range": "A1:B2", "alpha": "x", "groupedBy": 5}),
|
||
json!({"tool": "descriptive", "range": "A1:A2", "confidence": 250}),
|
||
json!({"tool": "anovaTwoFactorWithout", "range": "A1:B2"}),
|
||
json!({"tool": "anovaTwoFactorWithout", "range": "A1:A2"}),
|
||
json!({"tool": "anovaTwoFactorWith", "range": "A1:B2", "rowsPerSample": 2}),
|
||
json!({"tool": "anovaTwoFactorWith", "range": "A1:B2", "rowsPerSample": 1e300}),
|
||
json!({"tool": "zTest", "range1": "A1", "range2": "A2", "variance1": "x", "variance2": -1}),
|
||
] {
|
||
assert!(s.execute("data.analysisToolPak", p.clone()).is_err(), "{p}");
|
||
}
|
||
// Collinear X columns are an error, not a NaN report.
|
||
s.execute("range.setValues", json!({"range": "D1", "values": [[1, 1, 2], [2, 2, 4], [3, 3, 6], [5, 4, 8]]})).unwrap();
|
||
assert!(s.execute("data.analysisToolPak", json!({"tool": "regression", "inputY": "D1:D4", "inputX": "E1:F4"})).is_err());
|
||
// Constant data gives #DIV/0! cells, not a crash.
|
||
s.execute("range.setValues", json!({"range": "H1", "values": [[4], [4], [4]]})).unwrap();
|
||
s.execute("data.analysisToolPak", json!({"tool": "descriptive", "range": "H1:H3", "output": "J1"})).unwrap();
|
||
assert_eq!(v(&s, "K9"), Value::Error(gridcraft_core::CellError::Div0));
|
||
}
|
||
}
|