mirror of
https://github.com/storytold/gridcraft.git
synced 2026-10-27 20:26:48 -06:00
refactor(functions): move stat.rs items to stat/
This commit is contained in:
aiongg
+20
-1585
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,547 @@
|
||||
use super::*;
|
||||
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
// Basic helpers over number lists
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
|
||||
pub(super) fn mean(x: &[f64]) -> R<f64> {
|
||||
if x.is_empty() {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
Ok(x.iter().sum::<f64>() / x.len() as f64)
|
||||
}
|
||||
|
||||
pub(super) fn devsq(x: &[f64]) -> f64 {
|
||||
if x.is_empty() {
|
||||
return 0.0;
|
||||
}
|
||||
let m = x.iter().sum::<f64>() / x.len() as f64;
|
||||
x.iter().map(|v| (v - m) * (v - m)).sum()
|
||||
}
|
||||
|
||||
pub(super) fn var_s(x: &[f64]) -> R<f64> {
|
||||
if x.len() < 2 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
Ok(devsq(x) / (x.len() - 1) as f64)
|
||||
}
|
||||
|
||||
pub(super) fn var_p(x: &[f64]) -> R<f64> {
|
||||
if x.is_empty() {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
Ok(devsq(x) / x.len() as f64)
|
||||
}
|
||||
|
||||
fn sorted(mut v: Vec<f64>) -> Vec<f64> {
|
||||
v.sort_by(f64::total_cmp);
|
||||
v
|
||||
}
|
||||
|
||||
pub(super) fn median(v: Vec<f64>) -> R<f64> {
|
||||
if v.is_empty() {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let s = sorted(v);
|
||||
let n = s.len();
|
||||
let hi = s.get(n / 2).copied().unwrap_or(0.0);
|
||||
if n % 2 == 1 { Ok(hi) } else { Ok(s.get(n / 2 - 1).copied().unwrap_or(0.0).midpoint(hi)) }
|
||||
}
|
||||
|
||||
/// Values occurring at least twice, ordered by first occurrence, with their (equal, maximal) count.
|
||||
fn modes(v: &[f64]) -> Vec<f64> {
|
||||
let mut counts: Vec<(f64, usize)> = Vec::new();
|
||||
for &x in v {
|
||||
if let Some(e) = counts.iter_mut().find(|e| e.0 == x) {
|
||||
e.1 += 1;
|
||||
} else {
|
||||
counts.push((x, 1));
|
||||
}
|
||||
}
|
||||
let best = counts.iter().map(|e| e.1).max().unwrap_or(0);
|
||||
if best < 2 {
|
||||
return Vec::new();
|
||||
}
|
||||
counts.into_iter().filter(|e| e.1 == best).map(|e| e.0).collect()
|
||||
}
|
||||
|
||||
pub(super) fn mode_sngl(v: &[f64]) -> R<f64> {
|
||||
modes(v).first().copied().ok_or(CellError::NA)
|
||||
}
|
||||
|
||||
/// Numbers for MODE / MODE.SNGL / MODE.MULT. Inside multi-cell ranges and arrays only numbers
|
||||
/// count, like [`numbers`]. Scalar arguments differ from SUM: in Excel 365, text and logical
|
||||
/// values are skipped whether typed or referenced (even numeric text such as `"1"`, so
|
||||
/// `MODE(1,"1",2,3)` is `#N/A`, and `MODE(1,TRUE,TRUE,2)` is `#N/A`), but a single-cell reference
|
||||
/// to a blank cell is `#VALUE!` (`MODE(B9,2,D7)` with D7 blank). Arguments are processed left to
|
||||
/// right, so the first error wins. The results cached in POI's FormulaEvalTestData come from an
|
||||
/// older Excel that answered `#VALUE!` for typed text; they are not the reference here.
|
||||
pub(super) fn mode_numbers(args: &[Arg]) -> R<Vec<f64>> {
|
||||
let mut out = Vec::new();
|
||||
for a in args {
|
||||
match &a.value {
|
||||
Value::Array(arr) => {
|
||||
for v in arr.iter() {
|
||||
match v {
|
||||
Value::Number(n) => out.push(*n),
|
||||
Value::Error(e) => return Err(*e),
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
}
|
||||
Value::Number(n) => out.push(*n),
|
||||
Value::Error(e) => return Err(*e),
|
||||
Value::Text(_) | Value::Bool(_) => {}
|
||||
Value::Empty => return Err(CellError::Value),
|
||||
}
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
fn product(v: &[f64]) -> f64 {
|
||||
if v.is_empty() { 0.0 } else { v.iter().product() }
|
||||
}
|
||||
|
||||
pub(super) fn min_of(v: &[f64]) -> f64 {
|
||||
if v.is_empty() { 0.0 } else { v.iter().copied().fold(f64::INFINITY, f64::min) }
|
||||
}
|
||||
|
||||
pub(super) fn max_of(v: &[f64]) -> f64 {
|
||||
if v.is_empty() { 0.0 } else { v.iter().copied().fold(f64::NEG_INFINITY, f64::max) }
|
||||
}
|
||||
|
||||
/// k-th largest (k is rounded up like Excel).
|
||||
fn kth(v: Vec<f64>, k: f64, largest: bool) -> R<f64> {
|
||||
let k = k.ceil();
|
||||
if k.is_nan() || k < 1.0 || k > v.len() as f64 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let s = sorted(v);
|
||||
let k = k as usize;
|
||||
let idx = if largest { s.len() - k } else { k - 1 };
|
||||
s.get(idx).copied().ok_or(CellError::Num)
|
||||
}
|
||||
|
||||
fn interp(s: &[f64], pos: f64) -> R<f64> {
|
||||
let i = pos.floor();
|
||||
let frac = pos - i;
|
||||
let i = i as usize;
|
||||
let lo = *s.get(i).ok_or(CellError::Num)?;
|
||||
match s.get(i + 1) {
|
||||
Some(hi) if frac > 0.0 => Ok(lo + frac * (hi - lo)),
|
||||
_ => Ok(lo),
|
||||
}
|
||||
}
|
||||
|
||||
pub(super) fn percentile_inc(v: Vec<f64>, k: f64) -> R<f64> {
|
||||
if v.is_empty() || !(0.0..=1.0).contains(&k) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let s = sorted(v);
|
||||
interp(&s, k * (s.len() - 1) as f64)
|
||||
}
|
||||
|
||||
pub(super) fn percentile_exc(v: Vec<f64>, k: f64) -> R<f64> {
|
||||
let n = v.len() as f64;
|
||||
if v.is_empty() || !(k > 0.0 && k < 1.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let pos = k * (n + 1.0);
|
||||
if pos < 1.0 || pos > n {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let s = sorted(v);
|
||||
interp(&s, pos - 1.0)
|
||||
}
|
||||
|
||||
pub(super) fn quartile_inc(v: Vec<f64>, q: f64) -> R<f64> {
|
||||
let q = q.trunc();
|
||||
if !(0.0..=4.0).contains(&q) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
percentile_inc(v, q / 4.0)
|
||||
}
|
||||
|
||||
pub(super) fn quartile_exc(v: Vec<f64>, q: f64) -> R<f64> {
|
||||
let q = q.trunc();
|
||||
if !(1.0..=3.0).contains(&q) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
percentile_exc(v, q / 4.0)
|
||||
}
|
||||
|
||||
pub(super) fn arg_nums(a: &[Arg], i: usize) -> R<Vec<f64>> {
|
||||
array_numbers(&arg(a, i)?.value)
|
||||
}
|
||||
|
||||
/// Pairs of numbers at matching positions of two arrays (pairs with a non-number skipped).
|
||||
/// Different sizes give `#N/A`; errors propagate.
|
||||
fn pairs(x: &Value, y: &Value) -> R<(Vec<f64>, Vec<f64>)> {
|
||||
let xa = as_array(x);
|
||||
let ya = as_array(y);
|
||||
if xa.data.len() != ya.data.len() {
|
||||
return Err(CellError::NA);
|
||||
}
|
||||
let mut xs = Vec::new();
|
||||
let mut ys = Vec::new();
|
||||
for (a, b) in xa.data.iter().zip(ya.data.iter()) {
|
||||
if let Value::Error(e) = a {
|
||||
return Err(*e);
|
||||
}
|
||||
if let Value::Error(e) = b {
|
||||
return Err(*e);
|
||||
}
|
||||
if let (Value::Number(p), Value::Number(q)) = (a, b) {
|
||||
xs.push(*p);
|
||||
ys.push(*q);
|
||||
}
|
||||
}
|
||||
Ok((xs, ys))
|
||||
}
|
||||
|
||||
pub(super) fn arg_pairs(a: &[Arg], i: usize, j: usize) -> R<(Vec<f64>, Vec<f64>)> {
|
||||
pairs(&arg(a, i)?.value, &arg(a, j)?.value)
|
||||
}
|
||||
|
||||
/// (sum of cross deviations, devsq x, devsq y, mean x, mean y).
|
||||
pub(super) fn moments(x: &[f64], y: &[f64]) -> R<(f64, f64, f64, f64, f64)> {
|
||||
let mx = mean(x)?;
|
||||
let my = mean(y)?;
|
||||
let mut sxy = 0.0;
|
||||
let mut sxx = 0.0;
|
||||
let mut syy = 0.0;
|
||||
for (a, b) in x.iter().zip(y) {
|
||||
sxy += (a - mx) * (b - my);
|
||||
sxx += (a - mx) * (a - mx);
|
||||
syy += (b - my) * (b - my);
|
||||
}
|
||||
Ok((sxy, sxx, syy, mx, my))
|
||||
}
|
||||
|
||||
pub(super) fn correl(x: &[f64], y: &[f64]) -> R<f64> {
|
||||
if x.len() < 2 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let (sxy, sxx, syy, _, _) = moments(x, y)?;
|
||||
if sxx == 0.0 || syy == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
Ok(sxy / (sxx * syy).sqrt())
|
||||
}
|
||||
|
||||
pub(super) fn slope(y: &[f64], x: &[f64]) -> R<f64> {
|
||||
let (sxy, sxx, _, _, _) = moments(x, y)?;
|
||||
if sxx == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
Ok(sxy / sxx)
|
||||
}
|
||||
|
||||
pub(super) fn intercept(y: &[f64], x: &[f64]) -> R<f64> {
|
||||
let (sxy, sxx, _, mx, my) = moments(x, y)?;
|
||||
if sxx == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
Ok(my - sxy / sxx * mx)
|
||||
}
|
||||
|
||||
fn trunc_digits(r: f64, sig: f64) -> f64 {
|
||||
let f = 10f64.powi(sig as i32);
|
||||
if !f.is_finite() {
|
||||
return r;
|
||||
}
|
||||
((r * f) + 1e-9).floor() / f
|
||||
}
|
||||
|
||||
fn col_value(v: Vec<f64>) -> R<Value> {
|
||||
let n = v.len();
|
||||
array_val(n, 1, v.into_iter().map(Value::Number).collect())
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
// Aggregates
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
|
||||
pub(super) fn count_impl(a: &[Arg]) -> f64 {
|
||||
let mut c = 0usize;
|
||||
for x in a {
|
||||
match &x.value {
|
||||
Value::Array(arr) => c += arr.iter().filter(|v| v.is_number()).count(),
|
||||
Value::Number(_) => c += 1,
|
||||
_ if x.from_ref => {}
|
||||
Value::Bool(_) => c += 1,
|
||||
Value::Text(t) if gridcraft_core::parse::parse_number_text(t).is_some() => c += 1,
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
c as f64
|
||||
}
|
||||
|
||||
pub(super) fn counta_impl(a: &[Arg]) -> f64 {
|
||||
let mut c = 0usize;
|
||||
for x in a {
|
||||
match &x.value {
|
||||
Value::Array(arr) => c += arr.iter().filter(|v| !v.is_empty()).count(),
|
||||
Value::Empty if x.from_ref => {}
|
||||
_ => c += 1,
|
||||
}
|
||||
}
|
||||
c as f64
|
||||
}
|
||||
|
||||
pub(super) fn countblank(a: &[Arg], _: C) -> R<Value> {
|
||||
let v = &arg(a, 0)?.value;
|
||||
let blank = |x: &Value| matches!(x, Value::Empty) || matches!(x, Value::Text(t) if t.is_empty());
|
||||
let c = match v {
|
||||
Value::Array(arr) => arr.iter().filter(|x| blank(x)).count(),
|
||||
x => usize::from(blank(x)),
|
||||
};
|
||||
num_val(c as f64)
|
||||
}
|
||||
|
||||
pub(super) fn mode_mult(a: &[Arg], _: C) -> R<Value> {
|
||||
let v = mode_numbers(a)?;
|
||||
let m = modes(&v);
|
||||
if m.is_empty() {
|
||||
return Err(CellError::NA);
|
||||
}
|
||||
col_value(m)
|
||||
}
|
||||
|
||||
pub(super) fn large(a: &[Arg], _: C) -> R<Value> {
|
||||
num_val(kth(arg_nums(a, 0)?, num(a, 1)?, true)?)
|
||||
}
|
||||
|
||||
pub(super) fn small(a: &[Arg], _: C) -> R<Value> {
|
||||
num_val(kth(arg_nums(a, 0)?, num(a, 1)?, false)?)
|
||||
}
|
||||
|
||||
pub(super) fn rank_impl(a: &[Arg], avg: bool) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let v = arg_nums(a, 1)?;
|
||||
let asc = opt_num(a, 2, 0.0)? != 0.0;
|
||||
let eq = v.iter().filter(|&&y| y == x).count();
|
||||
if eq == 0 {
|
||||
return Err(CellError::NA);
|
||||
}
|
||||
let before = v.iter().filter(|&&y| if asc { y < x } else { y > x }).count();
|
||||
let r = before as f64 + 1.0;
|
||||
num_val(if avg { r + (eq as f64 - 1.0) / 2.0 } else { r })
|
||||
}
|
||||
|
||||
pub(super) fn percentrank_impl(a: &[Arg], exc: bool) -> R<Value> {
|
||||
let s = sorted(arg_nums(a, 0)?);
|
||||
let x = num(a, 1)?;
|
||||
let sig = opt_num(a, 2, 3.0)?.trunc();
|
||||
if sig < 1.0 || s.is_empty() {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let (first, last) = (s.first().copied().unwrap_or(0.0), s.last().copied().unwrap_or(0.0));
|
||||
if x < first || x > last {
|
||||
return Err(CellError::NA);
|
||||
}
|
||||
let n = s.len() as f64;
|
||||
let rank_of = |v: f64| -> f64 {
|
||||
let less = s.iter().filter(|&&y| y < v).count() as f64;
|
||||
if exc {
|
||||
(less + 1.0) / (n + 1.0)
|
||||
} else if n <= 1.0 {
|
||||
1.0
|
||||
} else {
|
||||
less / (n - 1.0)
|
||||
}
|
||||
};
|
||||
let r = if s.contains(&x) {
|
||||
rank_of(x)
|
||||
} else {
|
||||
let lo = s.iter().copied().filter(|&y| y < x).fold(f64::NEG_INFINITY, f64::max);
|
||||
let hi = s.iter().copied().filter(|&y| y > x).fold(f64::INFINITY, f64::min);
|
||||
let (rl, rh) = (rank_of(lo), rank_of(hi));
|
||||
rl + (x - lo) / (hi - lo) * (rh - rl)
|
||||
};
|
||||
num_val(trunc_digits(r, sig))
|
||||
}
|
||||
|
||||
pub(super) fn trimmean(a: &[Arg], _: C) -> R<Value> {
|
||||
let v = arg_nums(a, 0)?;
|
||||
let p = num(a, 1)?;
|
||||
if !(0.0..1.0).contains(&p) || v.is_empty() {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let n = v.len();
|
||||
let per_side = ((n as f64 * p) / 2.0).floor() as usize;
|
||||
let s = sorted(v);
|
||||
let kept = s.get(per_side..n - per_side).unwrap_or(&[]);
|
||||
num_val(mean(kept)?)
|
||||
}
|
||||
|
||||
pub(super) fn kurt(v: &[f64]) -> R<f64> {
|
||||
let n = v.len() as f64;
|
||||
if v.len() < 4 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let m = mean(v)?;
|
||||
let s = var_s(v)?.sqrt();
|
||||
if s == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let sum4: f64 = v.iter().map(|x| ((x - m) / s).powi(4)).sum();
|
||||
Ok(n * (n + 1.0) / ((n - 1.0) * (n - 2.0) * (n - 3.0)) * sum4 - 3.0 * (n - 1.0).powi(2) / ((n - 2.0) * (n - 3.0)))
|
||||
}
|
||||
|
||||
pub(super) fn skew(v: &[f64]) -> R<f64> {
|
||||
let n = v.len() as f64;
|
||||
if v.len() < 3 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let m = mean(v)?;
|
||||
let s = var_s(v)?.sqrt();
|
||||
if s == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let sum3: f64 = v.iter().map(|x| ((x - m) / s).powi(3)).sum();
|
||||
Ok(n / ((n - 1.0) * (n - 2.0)) * sum3)
|
||||
}
|
||||
|
||||
pub(super) fn skew_p(v: &[f64]) -> R<f64> {
|
||||
let m = mean(v)?;
|
||||
let s = var_p(v)?.sqrt();
|
||||
if s == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
Ok(v.iter().map(|x| ((x - m) / s).powi(3)).sum::<f64>() / v.len() as f64)
|
||||
}
|
||||
|
||||
pub(super) fn geomean(v: &[f64]) -> R<f64> {
|
||||
if v.is_empty() || v.iter().any(|&x| x <= 0.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
Ok((v.iter().map(|x| x.ln()).sum::<f64>() / v.len() as f64).exp())
|
||||
}
|
||||
|
||||
pub(super) fn harmean(v: &[f64]) -> R<f64> {
|
||||
if v.is_empty() || v.iter().any(|&x| x <= 0.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
Ok(v.len() as f64 / v.iter().map(|x| 1.0 / x).sum::<f64>())
|
||||
}
|
||||
|
||||
pub(super) fn avedev(v: &[f64]) -> R<f64> {
|
||||
if v.is_empty() {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let m = mean(v)?;
|
||||
Ok(v.iter().map(|x| (x - m).abs()).sum::<f64>() / v.len() as f64)
|
||||
}
|
||||
|
||||
pub(super) fn frequency(a: &[Arg], _: C) -> R<Value> {
|
||||
let data = arg_nums(a, 0)?;
|
||||
let bins = arg_nums(a, 1)?;
|
||||
if bins.is_empty() {
|
||||
return col_value(vec![data.len() as f64]);
|
||||
}
|
||||
let mut order: Vec<usize> = (0..bins.len()).collect();
|
||||
order.sort_by(|&i, &j| bins.get(i).copied().unwrap_or(0.0).total_cmp(&bins.get(j).copied().unwrap_or(0.0)));
|
||||
let mut counts = vec![0.0; bins.len() + 1];
|
||||
for x in data {
|
||||
let slot = order.iter().find(|&&i| x <= bins.get(i).copied().unwrap_or(f64::INFINITY)).copied().unwrap_or(bins.len());
|
||||
if let Some(c) = counts.get_mut(slot) {
|
||||
*c += 1.0;
|
||||
}
|
||||
}
|
||||
col_value(counts)
|
||||
}
|
||||
|
||||
pub(super) fn prob(a: &[Arg], _: C) -> R<Value> {
|
||||
let (xs, ps) = arg_pairs(a, 0, 1)?;
|
||||
let lo = num(a, 2)?;
|
||||
let hi = if has(a, 3) { num(a, 3)? } else { lo };
|
||||
if ps.iter().any(|&p| p <= 0.0 || p > 1.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let total: f64 = ps.iter().sum();
|
||||
if (total - 1.0).abs() > 1e-7 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(xs.iter().zip(&ps).filter(|(x, _)| **x >= lo && **x <= hi).map(|(_, p)| p).sum())
|
||||
}
|
||||
|
||||
pub(super) fn standardize(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
|
||||
if sd <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val((x - m) / sd)
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
// SUBTOTAL / AGGREGATE support
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
|
||||
/// Evaluates SUBTOTAL/AGGREGATE function `function_num` over already-filtered cell values.
|
||||
///
|
||||
/// Codes 1–11 (and 101–111) are AVERAGE, COUNT, COUNTA, MAX, MIN, PRODUCT, STDEV.S, STDEV.P,
|
||||
/// SUM, VAR.S, VAR.P; 12 = MEDIAN, 13 = MODE.SNGL. Codes 14–19 need `k`: use
|
||||
/// [`aggregate_values_k`] (they give `#VALUE!` here). Values that are not numbers are ignored
|
||||
/// (COUNTA counts non-empty values); the first error value present is returned.
|
||||
pub fn aggregate_values(function_num: u32, values: &[Value]) -> Value {
|
||||
crate::util::finish(aggregate_impl(function_num, values, None))
|
||||
}
|
||||
|
||||
/// AGGREGATE codes 14 = LARGE, 15 = SMALL, 16 = PERCENTILE.INC, 17 = QUARTILE.INC,
|
||||
/// 18 = PERCENTILE.EXC, 19 = QUARTILE.EXC with their `k` argument; codes 1–13 delegate to
|
||||
/// [`aggregate_values`].
|
||||
pub fn aggregate_values_k(function_num: u32, values: &[Value], k: &Value) -> Value {
|
||||
crate::util::finish(aggregate_impl(function_num, values, Some(k)))
|
||||
}
|
||||
|
||||
fn aggregate_impl(function_num: u32, values: &[Value], k: Option<&Value>) -> R<Value> {
|
||||
let code = if (101..=111).contains(&function_num) { function_num - 100 } else { function_num };
|
||||
let mut xs = Vec::new();
|
||||
let mut nonempty = 0usize;
|
||||
let mut push = |v: &Value| -> R<()> {
|
||||
match v {
|
||||
Value::Number(n) => xs.push(*n),
|
||||
Value::Error(e) => return Err(*e),
|
||||
_ => {}
|
||||
}
|
||||
if !v.is_empty() {
|
||||
nonempty += 1;
|
||||
}
|
||||
Ok(())
|
||||
};
|
||||
for v in values {
|
||||
match v {
|
||||
Value::Array(arr) => arr.iter().try_for_each(&mut push)?,
|
||||
other => push(other)?,
|
||||
}
|
||||
}
|
||||
let k = match (code, k) {
|
||||
(14..=19, Some(k)) => k.scalar().to_number()?,
|
||||
(14..=19, None) => return Err(CellError::Value),
|
||||
_ => 0.0,
|
||||
};
|
||||
let r = match code {
|
||||
1 => mean(&xs)?,
|
||||
2 => xs.len() as f64,
|
||||
3 => nonempty as f64,
|
||||
4 => max_of(&xs),
|
||||
5 => min_of(&xs),
|
||||
6 => product(&xs),
|
||||
7 => var_s(&xs)?.sqrt(),
|
||||
8 => var_p(&xs)?.sqrt(),
|
||||
9 => xs.iter().sum(),
|
||||
10 => var_s(&xs)?,
|
||||
11 => var_p(&xs)?,
|
||||
12 => median(xs)?,
|
||||
13 => mode_sngl(&xs)?,
|
||||
14 => kth(xs, k, true)?,
|
||||
15 => kth(xs, k, false)?,
|
||||
16 => percentile_inc(xs, k)?,
|
||||
17 => quartile_inc(xs, k)?,
|
||||
18 => percentile_exc(xs, k)?,
|
||||
19 => quartile_exc(xs, k)?,
|
||||
_ => return Err(CellError::Value),
|
||||
};
|
||||
num_val(r)
|
||||
}
|
||||
@@ -0,0 +1,702 @@
|
||||
use super::descriptive::{arg_nums, arg_pairs, mean, var_s};
|
||||
use super::*;
|
||||
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
// Distributions
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
|
||||
fn ln_comb(n: f64, k: f64) -> f64 {
|
||||
ln_gamma(n + 1.0) - ln_gamma(k + 1.0) - ln_gamma(n - k + 1.0)
|
||||
}
|
||||
|
||||
fn t_upper(t: f64, v: f64) -> f64 {
|
||||
// P(T > t) for t ≥ 0.
|
||||
0.5 * reg_inc_beta(v / (v + t * t), v / 2.0, 0.5)
|
||||
}
|
||||
|
||||
fn t_cdf(t: f64, v: f64) -> f64 {
|
||||
if t >= 0.0 { 1.0 - t_upper(t, v) } else { t_upper(-t, v) }
|
||||
}
|
||||
|
||||
fn t_pdf(t: f64, v: f64) -> f64 {
|
||||
(ln_gamma((v + 1.0) / 2.0) - ln_gamma(v / 2.0) - 0.5 * (v * std::f64::consts::PI).ln() - (v + 1.0) / 2.0 * (1.0 + t * t / v).ln()).exp()
|
||||
}
|
||||
|
||||
/// t ≥ 0 with P(T > t) = q, q in (0, 0.5].
|
||||
fn t_inv_upper(q: f64, v: f64) -> f64 {
|
||||
if q >= 0.5 {
|
||||
return 0.0;
|
||||
}
|
||||
invert(|t| t_upper(t, v), q, 0.0, 1.0, 1e60)
|
||||
}
|
||||
|
||||
fn chi_inv(p: f64, v: f64) -> f64 {
|
||||
if p <= 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
invert(|x| reg_gamma_p(v / 2.0, x / 2.0), p, 0.0, v.max(1.0), 1e60)
|
||||
}
|
||||
|
||||
fn chi_inv_rt(q: f64, v: f64) -> f64 {
|
||||
if q >= 1.0 {
|
||||
return 0.0;
|
||||
}
|
||||
invert(|x| reg_gamma_q(v / 2.0, x / 2.0), q, 0.0, v.max(1.0), 1e60)
|
||||
}
|
||||
|
||||
fn f_cdf(x: f64, d1: f64, d2: f64) -> f64 {
|
||||
if x <= 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
reg_inc_beta(d1 * x / (d1 * x + d2), d1 / 2.0, d2 / 2.0)
|
||||
}
|
||||
|
||||
fn f_upper(x: f64, d1: f64, d2: f64) -> f64 {
|
||||
if x <= 0.0 {
|
||||
return 1.0;
|
||||
}
|
||||
reg_inc_beta(d2 / (d2 + d1 * x), d2 / 2.0, d1 / 2.0)
|
||||
}
|
||||
|
||||
fn f_pdf(x: f64, d1: f64, d2: f64) -> f64 {
|
||||
if x == 0.0 {
|
||||
return if d1 < 2.0 {
|
||||
f64::INFINITY
|
||||
} else if d1 == 2.0 {
|
||||
1.0
|
||||
} else {
|
||||
0.0
|
||||
};
|
||||
}
|
||||
((d1 / 2.0) * (d1 / d2).ln() + (d1 / 2.0 - 1.0) * x.ln() - ((d1 + d2) / 2.0) * (1.0 + d1 * x / d2).ln() - ln_beta(d1 / 2.0, d2 / 2.0)).exp()
|
||||
}
|
||||
|
||||
fn binom_pmf(k: f64, n: f64, p: f64) -> f64 {
|
||||
if k < 0.0 || k > n {
|
||||
return 0.0;
|
||||
}
|
||||
if p == 0.0 {
|
||||
return if k == 0.0 { 1.0 } else { 0.0 };
|
||||
}
|
||||
if p == 1.0 {
|
||||
return if k == n { 1.0 } else { 0.0 };
|
||||
}
|
||||
(ln_comb(n, k) + k * p.ln() + (n - k) * (1.0 - p).ln()).exp()
|
||||
}
|
||||
|
||||
fn binom_cdf(k: f64, n: f64, p: f64) -> f64 {
|
||||
if k < 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
if k >= n {
|
||||
return 1.0;
|
||||
}
|
||||
if n <= 10_000.0 {
|
||||
let mut s = 0.0;
|
||||
let mut i = 0.0;
|
||||
while i <= k {
|
||||
s += binom_pmf(i, n, p);
|
||||
i += 1.0;
|
||||
}
|
||||
return s.min(1.0);
|
||||
}
|
||||
if p == 0.0 {
|
||||
return 1.0;
|
||||
}
|
||||
if p == 1.0 {
|
||||
return 0.0;
|
||||
}
|
||||
reg_inc_beta(1.0 - p, n - k, k + 1.0)
|
||||
}
|
||||
|
||||
fn poisson_pmf(x: f64, m: f64) -> f64 {
|
||||
if m == 0.0 {
|
||||
return if x == 0.0 { 1.0 } else { 0.0 };
|
||||
}
|
||||
(-m + x * m.ln() - ln_gamma(x + 1.0)).exp()
|
||||
}
|
||||
|
||||
fn poisson_cdf(x: f64, m: f64) -> f64 {
|
||||
if m == 0.0 {
|
||||
return 1.0;
|
||||
}
|
||||
if x <= 1000.0 {
|
||||
let mut s = 0.0;
|
||||
let mut i = 0.0;
|
||||
while i <= x {
|
||||
s += poisson_pmf(i, m);
|
||||
i += 1.0;
|
||||
}
|
||||
return s.min(1.0);
|
||||
}
|
||||
reg_gamma_q(x + 1.0, m)
|
||||
}
|
||||
|
||||
pub(super) fn norm_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, m, sd, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
|
||||
if sd <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let z = (x - m) / sd;
|
||||
num_val(if cum { norm_cdf(z) } else { norm_pdf(z) / sd })
|
||||
}
|
||||
|
||||
pub(super) fn norm_inv_fn(a: &[Arg], _: C) -> R<Value> {
|
||||
let (p, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
|
||||
if !(p > 0.0 && p < 1.0) || sd <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(m + sd * norm_inv(p))
|
||||
}
|
||||
|
||||
pub(super) fn norm_s_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let z = num(a, 0)?;
|
||||
let cum = if a.len() > 1 { boolean(a, 1)? } else { true };
|
||||
num_val(if cum { norm_cdf(z) } else { norm_pdf(z) })
|
||||
}
|
||||
|
||||
pub(super) fn norm_s_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let p = num(a, 0)?;
|
||||
if !(p > 0.0 && p < 1.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(norm_inv(p))
|
||||
}
|
||||
|
||||
fn df_arg(a: &[Arg], i: usize) -> R<f64> {
|
||||
let v = num(a, i)?.trunc();
|
||||
if v < 1.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
Ok(v)
|
||||
}
|
||||
|
||||
pub(super) fn t_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
let cum = boolean(a, 2)?;
|
||||
num_val(if cum { t_cdf(x, v) } else { t_pdf(x, v) })
|
||||
}
|
||||
|
||||
pub(super) fn t_dist_2t(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
if x < 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(2.0 * t_upper(x, v))
|
||||
}
|
||||
|
||||
pub(super) fn t_dist_rt(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
num_val(1.0 - t_cdf(x, v))
|
||||
}
|
||||
|
||||
pub(super) fn tdist_legacy(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
let tails = num(a, 2)?.trunc();
|
||||
if x < 0.0 || !(tails == 1.0 || tails == 2.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(tails * t_upper(x, v))
|
||||
}
|
||||
|
||||
pub(super) fn t_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let p = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
if !(p > 0.0 && p < 1.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(if p < 0.5 { -t_inv_upper(p, v) } else { t_inv_upper(1.0 - p, v) })
|
||||
}
|
||||
|
||||
pub(super) fn t_inv_2t(a: &[Arg], _: C) -> R<Value> {
|
||||
let p = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
if !(p > 0.0 && p <= 1.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(t_inv_upper(p / 2.0, v))
|
||||
}
|
||||
|
||||
pub(super) fn t_test(a: &[Arg], _: C) -> R<Value> {
|
||||
let tails = num(a, 2)?.trunc();
|
||||
let ty = num(a, 3)?.trunc();
|
||||
if !(tails == 1.0 || tails == 2.0) || !(1.0..=3.0).contains(&ty) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let (t, v) = if ty == 1.0 {
|
||||
let (x, y) = arg_pairs(a, 0, 1)?;
|
||||
let d: Vec<f64> = x.iter().zip(&y).map(|(p, q)| p - q).collect();
|
||||
let n = d.len() as f64;
|
||||
let sd = var_s(&d)?.sqrt();
|
||||
if sd == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
(mean(&d)? / (sd / n.sqrt()), n - 1.0)
|
||||
} else {
|
||||
let x = arg_nums(a, 0)?;
|
||||
let y = arg_nums(a, 1)?;
|
||||
let (n1, n2) = (x.len() as f64, y.len() as f64);
|
||||
let (v1, v2) = (var_s(&x)?, var_s(&y)?);
|
||||
let diff = mean(&x)? - mean(&y)?;
|
||||
if ty == 2.0 {
|
||||
let df = n1 + n2 - 2.0;
|
||||
let sp = ((n1 - 1.0) * v1 + (n2 - 1.0) * v2) / df;
|
||||
let se = (sp * (1.0 / n1 + 1.0 / n2)).sqrt();
|
||||
if se == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
(diff / se, df)
|
||||
} else {
|
||||
let (a1, a2) = (v1 / n1, v2 / n2);
|
||||
let se = (a1 + a2).sqrt();
|
||||
if se == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let df = (a1 + a2).powi(2) / (a1 * a1 / (n1 - 1.0) + a2 * a2 / (n2 - 1.0));
|
||||
(diff / se, df)
|
||||
}
|
||||
};
|
||||
num_val(tails * t_upper(t.abs(), v))
|
||||
}
|
||||
|
||||
pub(super) fn chisq_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
let cum = boolean(a, 2)?;
|
||||
if x < 0.0 || v > 1e10 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if cum {
|
||||
return num_val(reg_gamma_p(v / 2.0, x / 2.0));
|
||||
}
|
||||
if x == 0.0 {
|
||||
return if v == 2.0 {
|
||||
num_val(0.5)
|
||||
} else if v < 2.0 {
|
||||
Err(CellError::Num)
|
||||
} else {
|
||||
num_val(0.0)
|
||||
};
|
||||
}
|
||||
num_val(((v / 2.0 - 1.0) * x.ln() - x / 2.0 - (v / 2.0) * 2f64.ln() - ln_gamma(v / 2.0)).exp())
|
||||
}
|
||||
|
||||
pub(super) fn chisq_dist_rt(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
if x < 0.0 || v > 1e10 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(reg_gamma_q(v / 2.0, x / 2.0))
|
||||
}
|
||||
|
||||
pub(super) fn chisq_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let p = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
if !(0.0..1.0).contains(&p) || v > 1e10 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(chi_inv(p, v))
|
||||
}
|
||||
|
||||
pub(super) fn chisq_inv_rt(a: &[Arg], _: C) -> R<Value> {
|
||||
let p = num(a, 0)?;
|
||||
let v = df_arg(a, 1)?;
|
||||
if !(p > 0.0 && p <= 1.0) || v > 1e10 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(chi_inv_rt(p, v))
|
||||
}
|
||||
|
||||
pub(super) fn chisq_test(a: &[Arg], _: C) -> R<Value> {
|
||||
let act = as_array(&arg(a, 0)?.value).into_owned();
|
||||
let exp = as_array(&arg(a, 1)?.value).into_owned();
|
||||
if act.rows != exp.rows || act.cols != exp.cols {
|
||||
return Err(CellError::NA);
|
||||
}
|
||||
let mut x2 = 0.0;
|
||||
for (o, e) in act.data.iter().zip(&exp.data) {
|
||||
if let Value::Error(err) = o {
|
||||
return Err(*err);
|
||||
}
|
||||
if let Value::Error(err) = e {
|
||||
return Err(*err);
|
||||
}
|
||||
if let (Value::Number(o), Value::Number(e)) = (o, e) {
|
||||
if *e == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
x2 += (o - e).powi(2) / e;
|
||||
}
|
||||
}
|
||||
let df = if act.rows > 1 && act.cols > 1 { ((act.rows - 1) * (act.cols - 1)) as f64 } else { (act.data.len().saturating_sub(1)) as f64 };
|
||||
if df < 1.0 {
|
||||
return Err(CellError::NA);
|
||||
}
|
||||
num_val(reg_gamma_q(df / 2.0, x2 / 2.0))
|
||||
}
|
||||
|
||||
fn f_args(a: &[Arg]) -> R<(f64, f64, f64)> {
|
||||
let x = num(a, 0)?;
|
||||
let d1 = num(a, 1)?.trunc();
|
||||
let d2 = num(a, 2)?.trunc();
|
||||
if d1 < 1.0 || d2 < 1.0 || d1 >= 1e10 || d2 >= 1e10 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
Ok((x, d1, d2))
|
||||
}
|
||||
|
||||
pub(super) fn f_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, d1, d2) = f_args(a)?;
|
||||
let cum = boolean(a, 3)?;
|
||||
if x < 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(if cum { f_cdf(x, d1, d2) } else { f_pdf(x, d1, d2) })
|
||||
}
|
||||
|
||||
pub(super) fn f_dist_rt(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, d1, d2) = f_args(a)?;
|
||||
if x < 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(f_upper(x, d1, d2))
|
||||
}
|
||||
|
||||
pub(super) fn f_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let (p, d1, d2) = f_args(a)?;
|
||||
if !(0.0..1.0).contains(&p) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if p == 0.0 {
|
||||
return num_val(0.0);
|
||||
}
|
||||
num_val(invert(|x| f_cdf(x, d1, d2), p, 0.0, 1.0, 1e60))
|
||||
}
|
||||
|
||||
pub(super) fn f_inv_rt(a: &[Arg], _: C) -> R<Value> {
|
||||
let (p, d1, d2) = f_args(a)?;
|
||||
if !(p > 0.0 && p <= 1.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if p == 1.0 {
|
||||
return num_val(0.0);
|
||||
}
|
||||
num_val(invert(|x| f_upper(x, d1, d2), p, 0.0, 1.0, 1e60))
|
||||
}
|
||||
|
||||
pub(super) fn f_test(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = arg_nums(a, 0)?;
|
||||
let y = arg_nums(a, 1)?;
|
||||
let (v1, v2) = (var_s(&x)?, var_s(&y)?);
|
||||
if v1 == 0.0 || v2 == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let f = v1 / v2;
|
||||
let (d1, d2) = ((x.len() - 1) as f64, (y.len() - 1) as f64);
|
||||
let c = f_cdf(f, d1, d2);
|
||||
num_val((2.0 * c.min(1.0 - c)).min(1.0))
|
||||
}
|
||||
|
||||
pub(super) fn binom_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let k = num(a, 0)?.trunc();
|
||||
let n = num(a, 1)?.trunc();
|
||||
let p = num(a, 2)?;
|
||||
let cum = boolean(a, 3)?;
|
||||
if k < 0.0 || k > n || !(0.0..=1.0).contains(&p) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(if cum { binom_cdf(k, n, p) } else { binom_pmf(k, n, p) })
|
||||
}
|
||||
|
||||
pub(super) fn binom_dist_range(a: &[Arg], _: C) -> R<Value> {
|
||||
let n = num(a, 0)?.trunc();
|
||||
let p = num(a, 1)?;
|
||||
let s1 = num(a, 2)?.trunc();
|
||||
let s2 = if has(a, 3) { num(a, 3)?.trunc() } else { s1 };
|
||||
if !(0.0..=1e15).contains(&n) || !(0.0..=1.0).contains(&p) || s1 < 0.0 || s1 > n || s2 < s1 || s2 > n {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if s2 - s1 > 10_000.0 {
|
||||
return num_val(binom_cdf(s2, n, p) - binom_cdf(s1 - 1.0, n, p));
|
||||
}
|
||||
let mut s = 0.0;
|
||||
let mut k = s1;
|
||||
while k <= s2 {
|
||||
s += binom_pmf(k, n, p);
|
||||
k += 1.0;
|
||||
}
|
||||
num_val(s)
|
||||
}
|
||||
|
||||
pub(super) fn binom_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let n = num(a, 0)?.trunc();
|
||||
let p = num(a, 1)?;
|
||||
let alpha = num(a, 2)?;
|
||||
if n < 0.0 || !(0.0..=1.0).contains(&p) || !(0.0..=1.0).contains(&alpha) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
// Smallest k with cdf(k) ≥ alpha, by bisection on the integer k.
|
||||
let (mut lo, mut hi) = (0.0f64, n);
|
||||
if binom_cdf(0.0, n, p) >= alpha {
|
||||
return num_val(0.0);
|
||||
}
|
||||
while hi - lo > 1.0 {
|
||||
let mid = ((lo + hi) / 2.0).floor();
|
||||
// Above 2^53 adjacent f64 values are more than 1 apart, so the gap can stop shrinking.
|
||||
if mid <= lo || mid >= hi {
|
||||
break;
|
||||
}
|
||||
if binom_cdf(mid, n, p) >= alpha {
|
||||
hi = mid;
|
||||
} else {
|
||||
lo = mid;
|
||||
}
|
||||
}
|
||||
num_val(hi)
|
||||
}
|
||||
|
||||
pub(super) fn poisson_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?.trunc();
|
||||
let m = num(a, 1)?;
|
||||
let cum = boolean(a, 2)?;
|
||||
if x < 0.0 || m < 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(if cum { poisson_cdf(x, m) } else { poisson_pmf(x, m) })
|
||||
}
|
||||
|
||||
pub(super) fn expon_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, l, cum) = (num(a, 0)?, num(a, 1)?, boolean(a, 2)?);
|
||||
if x < 0.0 || l <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(if cum { -(-l * x).exp_m1() } else { l * (-l * x).exp() })
|
||||
}
|
||||
|
||||
pub(super) fn gamma_fn(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let g = gamma(x);
|
||||
if !g.is_finite() {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(g)
|
||||
}
|
||||
|
||||
pub(super) fn gamma_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
|
||||
if x < 0.0 || al <= 0.0 || be <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if cum {
|
||||
return num_val(reg_gamma_p(al, x / be));
|
||||
}
|
||||
if x == 0.0 {
|
||||
return if al < 1.0 {
|
||||
Err(CellError::Num)
|
||||
} else if al == 1.0 {
|
||||
num_val(1.0 / be)
|
||||
} else {
|
||||
num_val(0.0)
|
||||
};
|
||||
}
|
||||
num_val(((al - 1.0) * x.ln() - x / be - al * be.ln() - ln_gamma(al)).exp())
|
||||
}
|
||||
|
||||
pub(super) fn gamma_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let (p, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
|
||||
if !(0.0..1.0).contains(&p) || al <= 0.0 || be <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if p == 0.0 {
|
||||
return num_val(0.0);
|
||||
}
|
||||
num_val(invert(|x| reg_gamma_p(al, x / be), p, 0.0, (al * be).max(1e-3), 1e60))
|
||||
}
|
||||
|
||||
pub(super) fn gammaln(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
if x <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(ln_gamma(x))
|
||||
}
|
||||
|
||||
fn beta_bounds(a: &[Arg], i: usize) -> R<(f64, f64)> {
|
||||
let lo = if has(a, i) { num(a, i)? } else { 0.0 };
|
||||
let hi = if has(a, i + 1) { num(a, i + 1)? } else { 1.0 };
|
||||
if lo >= hi {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
Ok((lo, hi))
|
||||
}
|
||||
|
||||
pub(super) fn beta_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
|
||||
let (lo, hi) = beta_bounds(a, 4)?;
|
||||
if al <= 0.0 || be <= 0.0 || x < lo || x > hi {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let z = (x - lo) / (hi - lo);
|
||||
if cum {
|
||||
return num_val(reg_inc_beta(z, al, be));
|
||||
}
|
||||
if (z == 0.0 && al < 1.0) || (z == 1.0 && be < 1.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let lp = if z == 0.0 || z == 1.0 {
|
||||
let base = if z == 0.0 { al } else { be };
|
||||
if base == 1.0 { -ln_beta(al, be) } else { f64::NEG_INFINITY }
|
||||
} else {
|
||||
(al - 1.0) * z.ln() + (be - 1.0) * (1.0 - z).ln() - ln_beta(al, be)
|
||||
};
|
||||
num_val(lp.exp() / (hi - lo))
|
||||
}
|
||||
|
||||
pub(super) fn betadist_legacy(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
|
||||
let (lo, hi) = beta_bounds(a, 3)?;
|
||||
if al <= 0.0 || be <= 0.0 || x < lo || x > hi {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(reg_inc_beta((x - lo) / (hi - lo), al, be))
|
||||
}
|
||||
|
||||
pub(super) fn beta_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let (p, al, be) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
|
||||
let (lo, hi) = beta_bounds(a, 3)?;
|
||||
if al <= 0.0 || be <= 0.0 || !(0.0..=1.0).contains(&p) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let z = if p == 0.0 {
|
||||
0.0
|
||||
} else if p == 1.0 {
|
||||
1.0
|
||||
} else {
|
||||
invert(|x| reg_inc_beta(x, al, be), p, 0.0, 1.0, 1.0)
|
||||
};
|
||||
num_val(lo + z * (hi - lo))
|
||||
}
|
||||
|
||||
pub(super) fn lognorm_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
|
||||
let cum = if a.len() > 3 { boolean(a, 3)? } else { true };
|
||||
if x <= 0.0 || sd <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let z = (x.ln() - m) / sd;
|
||||
num_val(if cum { norm_cdf(z) } else { norm_pdf(z) / (x * sd) })
|
||||
}
|
||||
|
||||
pub(super) fn lognorm_inv(a: &[Arg], _: C) -> R<Value> {
|
||||
let (p, m, sd) = (num(a, 0)?, num(a, 1)?, num(a, 2)?);
|
||||
if !(p > 0.0 && p < 1.0) || sd <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val((m + sd * norm_inv(p)).exp())
|
||||
}
|
||||
|
||||
pub(super) fn weibull_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let (x, al, be, cum) = (num(a, 0)?, num(a, 1)?, num(a, 2)?, boolean(a, 3)?);
|
||||
if x < 0.0 || al <= 0.0 || be <= 0.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let t = (x / be).powf(al);
|
||||
num_val(if cum { -(-t).exp_m1() } else { al / be.powf(al) * x.powf(al - 1.0) * (-t).exp() })
|
||||
}
|
||||
|
||||
pub(super) fn hypgeom_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let k = num(a, 0)?.trunc();
|
||||
let n = num(a, 1)?.trunc();
|
||||
let kk = num(a, 2)?.trunc();
|
||||
let nn = num(a, 3)?.trunc();
|
||||
let cum = if a.len() > 4 { boolean(a, 4)? } else { false };
|
||||
if k < 0.0 || n <= 0.0 || kk <= 0.0 || nn <= 0.0 || n > nn || kk > nn || k > n.min(kk) || k < (n - nn + kk).max(0.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
let pmf = |i: f64| (ln_comb(kk, i) + ln_comb(nn - kk, n - i) - ln_comb(nn, n)).exp();
|
||||
if !cum {
|
||||
return num_val(pmf(k));
|
||||
}
|
||||
let mut s = 0.0;
|
||||
let mut i = (n - nn + kk).max(0.0);
|
||||
if k - i > 1e6 || i + 1.0 == i {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
while i <= k {
|
||||
s += pmf(i);
|
||||
i += 1.0;
|
||||
}
|
||||
num_val(s.min(1.0))
|
||||
}
|
||||
|
||||
pub(super) fn negbinom_dist(a: &[Arg], _: C) -> R<Value> {
|
||||
let f = num(a, 0)?.trunc();
|
||||
let s = num(a, 1)?.trunc();
|
||||
let p = num(a, 2)?;
|
||||
let cum = if a.len() > 3 { boolean(a, 3)? } else { false };
|
||||
if f < 0.0 || s < 1.0 || !(0.0..=1.0).contains(&p) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if cum {
|
||||
if p == 0.0 {
|
||||
return num_val(0.0);
|
||||
}
|
||||
return num_val(reg_inc_beta(p, s, f + 1.0));
|
||||
}
|
||||
if p == 0.0 || p == 1.0 {
|
||||
return num_val(if p == 1.0 && f == 0.0 { 1.0 } else { 0.0 });
|
||||
}
|
||||
num_val((ln_comb(f + s - 1.0, s - 1.0) + s * p.ln() + f * (1.0 - p).ln()).exp())
|
||||
}
|
||||
|
||||
pub(super) fn confidence_norm(a: &[Arg], _: C) -> R<Value> {
|
||||
let (al, sd, n) = (num(a, 0)?, num(a, 1)?, num(a, 2)?.trunc());
|
||||
if !(al > 0.0 && al < 1.0) || sd <= 0.0 || n < 1.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(norm_inv(1.0 - al / 2.0) * sd / n.sqrt())
|
||||
}
|
||||
|
||||
pub(super) fn confidence_t(a: &[Arg], _: C) -> R<Value> {
|
||||
let (al, sd, n) = (num(a, 0)?, num(a, 1)?, num(a, 2)?.trunc());
|
||||
if !(al > 0.0 && al < 1.0) || sd <= 0.0 || n < 1.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
if n == 1.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
num_val(t_inv_upper(al / 2.0, n - 1.0) * sd / n.sqrt())
|
||||
}
|
||||
|
||||
pub(super) fn z_test(a: &[Arg], _: C) -> R<Value> {
|
||||
let v = arg_nums(a, 0)?;
|
||||
let x = num(a, 1)?;
|
||||
if v.is_empty() {
|
||||
return Err(CellError::NA);
|
||||
}
|
||||
let sd = if has(a, 2) { num(a, 2)? } else { var_s(&v)?.sqrt() };
|
||||
if sd <= 0.0 {
|
||||
return Err(if has(a, 2) { CellError::Num } else { CellError::Div0 });
|
||||
}
|
||||
let z = (mean(&v)? - x) / (sd / (v.len() as f64).sqrt());
|
||||
num_val(1.0 - norm_cdf(z))
|
||||
}
|
||||
|
||||
pub(super) fn fisher(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
if x <= -1.0 || x >= 1.0 {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
num_val(0.5 * ((1.0 + x) / (1.0 - x)).ln())
|
||||
}
|
||||
|
||||
pub(super) fn gauss(a: &[Arg], _: C) -> R<Value> {
|
||||
let z = num(a, 0)?;
|
||||
// Φ(z) − 0.5 = erf(z/√2)/2, accurate near zero.
|
||||
num_val(0.5 * erf(z / std::f64::consts::SQRT_2))
|
||||
}
|
||||
@@ -0,0 +1,342 @@
|
||||
use super::descriptive::{arg_pairs, intercept, moments, slope};
|
||||
use super::*;
|
||||
|
||||
pub(super) fn forecast(a: &[Arg], _: C) -> R<Value> {
|
||||
let x = num(a, 0)?;
|
||||
let (ys, xs) = arg_pairs(a, 1, 2)?;
|
||||
if ys.is_empty() {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
num_val(intercept(&ys, &xs)? + slope(&ys, &xs)? * x)
|
||||
}
|
||||
|
||||
pub(super) fn steyx(a: &[Arg], _: C) -> R<Value> {
|
||||
let (ys, xs) = arg_pairs(a, 0, 1)?;
|
||||
let n = ys.len() as f64;
|
||||
if ys.len() < 3 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let (sxy, sxx, syy, _, _) = moments(&xs, &ys)?;
|
||||
if sxx == 0.0 {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
num_val(((syy - sxy * sxy / sxx) / (n - 2.0)).max(0.0).sqrt())
|
||||
}
|
||||
|
||||
pub(super) fn covariance(a: &[Arg], sample: bool) -> R<Value> {
|
||||
let (xs, ys) = arg_pairs(a, 0, 1)?;
|
||||
let n = xs.len();
|
||||
if n == 0 || (sample && n < 2) {
|
||||
return Err(CellError::Div0);
|
||||
}
|
||||
let (sxy, ..) = moments(&xs, &ys)?;
|
||||
num_val(sxy / if sample { (n - 1) as f64 } else { n as f64 })
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
// Regression (LINEST, LOGEST, TREND, GROWTH)
|
||||
// ---------------------------------------------------------------------------------------------
|
||||
|
||||
struct Xy {
|
||||
y: Vec<f64>,
|
||||
/// n rows of k predictor values.
|
||||
x: Vec<Vec<f64>>,
|
||||
k: usize,
|
||||
y_is_row: bool,
|
||||
y_rows: usize,
|
||||
y_cols: usize,
|
||||
}
|
||||
|
||||
fn strict_num(v: &Value) -> R<f64> {
|
||||
match v {
|
||||
Value::Number(n) => Ok(*n),
|
||||
Value::Error(e) => Err(*e),
|
||||
_ => Err(CellError::Value),
|
||||
}
|
||||
}
|
||||
|
||||
/// Most predictor columns LINEST, LOGEST, TREND and GROWTH accept. The normal equations are a
|
||||
/// k×k matrix inverted in O(k³), so an unbounded k (one `SEQUENCE(1,16384)` argument) allocated
|
||||
/// gigabytes and could abort the process.
|
||||
pub(super) const MAX_REGRESSORS: usize = 255;
|
||||
|
||||
/// Most multiply-adds (rows × k²) a regression may spend building its normal equations.
|
||||
const MAX_REGRESSION_WORK: usize = 1_000_000_000;
|
||||
|
||||
/// Rejects regressions too large to fit: more than [`MAX_REGRESSORS`] predictors is `#VALUE!`,
|
||||
/// more than [`MAX_REGRESSION_WORK`] multiply-adds is `#NUM!`.
|
||||
pub(super) fn check_regression_size(n: usize, k: usize) -> R<()> {
|
||||
if k > MAX_REGRESSORS {
|
||||
return Err(CellError::Value);
|
||||
}
|
||||
match n.checked_mul(k).and_then(|nk| nk.checked_mul(k)) {
|
||||
Some(work) if work <= MAX_REGRESSION_WORK => Ok(()),
|
||||
_ => Err(CellError::Num),
|
||||
}
|
||||
}
|
||||
|
||||
fn parse_xy(y: &Value, x: Option<&Value>) -> R<Xy> {
|
||||
let ya = as_array(y);
|
||||
let n = ya.data.len();
|
||||
let y_is_row = ya.rows == 1 && ya.cols > 1;
|
||||
let yv: Vec<f64> = ya.data.iter().map(strict_num).collect::<R<_>>()?;
|
||||
let (xs, k) = match x {
|
||||
None => ((1..=n).map(|i| vec![i as f64]).collect(), 1),
|
||||
Some(xv) => {
|
||||
let xa = as_array(xv);
|
||||
if xa.rows == ya.rows && xa.cols == ya.cols {
|
||||
(xa.data.iter().map(|v| strict_num(v).map(|f| vec![f])).collect::<R<Vec<_>>>()?, 1)
|
||||
} else if !y_is_row && xa.rows == n {
|
||||
check_regression_size(n, xa.cols)?;
|
||||
let mut rows = Vec::with_capacity(n);
|
||||
for r in 0..n {
|
||||
rows.push((0..xa.cols).map(|c| strict_num(xa.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?);
|
||||
}
|
||||
(rows, xa.cols)
|
||||
} else if y_is_row && xa.cols == n {
|
||||
check_regression_size(n, xa.rows)?;
|
||||
let mut rows = Vec::with_capacity(n);
|
||||
for c in 0..n {
|
||||
rows.push((0..xa.rows).map(|r| strict_num(xa.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?);
|
||||
}
|
||||
(rows, xa.rows)
|
||||
} else {
|
||||
return Err(CellError::Ref);
|
||||
}
|
||||
}
|
||||
};
|
||||
Ok(Xy { y: yv, x: xs, k, y_is_row, y_rows: ya.rows, y_cols: ya.cols })
|
||||
}
|
||||
|
||||
struct Fit {
|
||||
coef: Vec<f64>,
|
||||
b: f64,
|
||||
se: Vec<f64>,
|
||||
se_b: Option<f64>,
|
||||
r2: R<f64>,
|
||||
sey: R<f64>,
|
||||
f: R<f64>,
|
||||
df: f64,
|
||||
ssreg: f64,
|
||||
ssresid: f64,
|
||||
}
|
||||
|
||||
/// Inverts a symmetric positive semi-definite matrix by Gauss–Jordan on the diagonal, dropping
|
||||
/// (zeroing) collinear variables. Returns the inverse and the dropped flags.
|
||||
fn sym_inverse(m: &[Vec<f64>]) -> (Vec<Vec<f64>>, Vec<bool>) {
|
||||
let p = m.len();
|
||||
let mut a: Vec<Vec<f64>> = m.to_vec();
|
||||
let mut inv: Vec<Vec<f64>> = (0..p).map(|i| (0..p).map(|j| if i == j { 1.0 } else { 0.0 }).collect()).collect();
|
||||
let mut dropped = vec![false; p];
|
||||
let diag: Vec<f64> = (0..p).map(|i| m.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0).abs()).collect();
|
||||
for j in 0..p {
|
||||
let piv = a.get(j).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
|
||||
let scale = diag.get(j).copied().unwrap_or(0.0);
|
||||
if piv.abs() <= 1e-12 * scale || piv.abs() < 1e-300 {
|
||||
if let Some(d) = dropped.get_mut(j) {
|
||||
*d = true;
|
||||
}
|
||||
for row in a.iter_mut() {
|
||||
if let Some(v) = row.get_mut(j) {
|
||||
*v = 0.0;
|
||||
}
|
||||
}
|
||||
if let Some(r) = a.get_mut(j) {
|
||||
r.iter_mut().for_each(|v| *v = 0.0);
|
||||
}
|
||||
if let Some(r) = inv.get_mut(j) {
|
||||
r.iter_mut().for_each(|v| *v = 0.0);
|
||||
}
|
||||
continue;
|
||||
}
|
||||
let arow: Vec<f64> = a.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default();
|
||||
let irow: Vec<f64> = inv.get(j).map(|r| r.iter().map(|v| v / piv).collect()).unwrap_or_default();
|
||||
for i in 0..p {
|
||||
if i == j {
|
||||
continue;
|
||||
}
|
||||
let f = a.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
|
||||
if f == 0.0 {
|
||||
continue;
|
||||
}
|
||||
if let Some(r) = a.get_mut(i) {
|
||||
r.iter_mut().zip(&arow).for_each(|(v, w)| *v -= f * w);
|
||||
}
|
||||
if let Some(r) = inv.get_mut(i) {
|
||||
r.iter_mut().zip(&irow).for_each(|(v, w)| *v -= f * w);
|
||||
}
|
||||
}
|
||||
if let Some(r) = a.get_mut(j) {
|
||||
*r = arow;
|
||||
}
|
||||
if let Some(r) = inv.get_mut(j) {
|
||||
*r = irow;
|
||||
}
|
||||
}
|
||||
(inv, dropped)
|
||||
}
|
||||
|
||||
fn regress(y: &[f64], x: &[Vec<f64>], k: usize, constant: bool) -> R<Fit> {
|
||||
let n = y.len();
|
||||
if n == 0 || x.len() != n {
|
||||
return Err(CellError::Value);
|
||||
}
|
||||
check_regression_size(n, k)?;
|
||||
let nf = n as f64;
|
||||
let (xm, ym) = if constant {
|
||||
let xm: Vec<f64> = (0..k).map(|j| x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)).sum::<f64>() / nf).collect();
|
||||
(xm, y.iter().sum::<f64>() / nf)
|
||||
} else {
|
||||
(vec![0.0; k], 0.0)
|
||||
};
|
||||
// Normal equations on (centred) data.
|
||||
let mut xtx = vec![vec![0.0; k]; k];
|
||||
let mut xty = vec![0.0; k];
|
||||
for (row, yv) in x.iter().zip(y) {
|
||||
let d: Vec<f64> = (0..k).map(|j| row.get(j).copied().unwrap_or(0.0) - xm.get(j).copied().unwrap_or(0.0)).collect();
|
||||
let yd = yv - ym;
|
||||
for (i, di) in d.iter().enumerate() {
|
||||
if let Some(t) = xty.get_mut(i) {
|
||||
*t += di * yd;
|
||||
}
|
||||
if let Some(r) = xtx.get_mut(i) {
|
||||
for (j, dj) in d.iter().enumerate() {
|
||||
if let Some(v) = r.get_mut(j) {
|
||||
*v += di * dj;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let (inv, dropped) = sym_inverse(&xtx);
|
||||
let coef: Vec<f64> = (0..k).map(|i| inv.get(i).map(|r| r.iter().zip(&xty).map(|(a, b)| a * b).sum()).unwrap_or(0.0)).collect();
|
||||
let b = if constant { ym - coef.iter().zip(&xm).map(|(c, m)| c * m).sum::<f64>() } else { 0.0 };
|
||||
let mut ssresid = 0.0;
|
||||
let mut sstot = 0.0;
|
||||
for (row, yv) in x.iter().zip(y) {
|
||||
let pred = b + coef.iter().zip(row).map(|(c, v)| c * v).sum::<f64>();
|
||||
ssresid += (yv - pred).powi(2);
|
||||
sstot += if constant { (yv - ym).powi(2) } else { yv * yv };
|
||||
}
|
||||
let used = dropped.iter().filter(|d| !**d).count() + usize::from(constant);
|
||||
let df = nf - used as f64;
|
||||
let ssreg = sstot - ssresid;
|
||||
let r2 = if sstot == 0.0 { Err(CellError::Num) } else { Ok(ssreg / sstot) };
|
||||
let sey2 = if df > 0.0 { Ok(ssresid / df) } else { Err(CellError::Num) };
|
||||
let sey = sey2.map(f64::sqrt);
|
||||
let v1 = if constant { nf - df - 1.0 } else { nf - df };
|
||||
let f = match sey2 {
|
||||
Ok(s2) if s2 > 0.0 && v1 > 0.0 => Ok((ssreg / v1) / s2),
|
||||
_ => Err(CellError::Num),
|
||||
};
|
||||
let s2 = sey2.unwrap_or(f64::NAN);
|
||||
let se: Vec<f64> = (0..k)
|
||||
.map(|i| if dropped.get(i).copied().unwrap_or(true) { 0.0 } else { (inv.get(i).and_then(|r| r.get(i)).copied().unwrap_or(0.0) * s2).sqrt() })
|
||||
.collect();
|
||||
let se_b = if constant {
|
||||
let mut q = 0.0;
|
||||
for (i, mi) in xm.iter().enumerate() {
|
||||
for (j, mj) in xm.iter().enumerate() {
|
||||
q += mi * mj * inv.get(i).and_then(|r| r.get(j)).copied().unwrap_or(0.0);
|
||||
}
|
||||
}
|
||||
Some((s2 * (1.0 / nf + q)).sqrt())
|
||||
} else {
|
||||
None
|
||||
};
|
||||
Ok(Fit { coef, b, se, se_b, r2, sey, f, df, ssreg, ssresid })
|
||||
}
|
||||
|
||||
fn const_arg(a: &[Arg], i: usize) -> R<bool> {
|
||||
if has(a, i) { boolean(a, i) } else { Ok(true) }
|
||||
}
|
||||
|
||||
fn opt_value(a: &[Arg], i: usize) -> Option<&Value> {
|
||||
if has(a, i) { a.get(i).map(|x| &x.value) } else { None }
|
||||
}
|
||||
|
||||
pub(super) fn linest_impl(a: &[Arg], log: bool) -> R<Value> {
|
||||
let mut xy = parse_xy(&arg(a, 0)?.value, opt_value(a, 1))?;
|
||||
if log {
|
||||
if xy.y.iter().any(|&v| v <= 0.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
xy.y.iter_mut().for_each(|v| *v = v.ln());
|
||||
}
|
||||
let constant = const_arg(a, 2)?;
|
||||
let stats = if has(a, 3) { boolean(a, 3)? } else { false };
|
||||
let fit = regress(&xy.y, &xy.x, xy.k, constant)?;
|
||||
let tr = |v: f64| if log { v.exp() } else { v };
|
||||
let cols = xy.k + 1;
|
||||
let mut data: Vec<Value> = Vec::new();
|
||||
for c in fit.coef.iter().rev() {
|
||||
data.push(Value::number(tr(*c)));
|
||||
}
|
||||
data.push(Value::number(tr(fit.b)));
|
||||
if stats {
|
||||
let rv = |r: R<f64>| r.map(Value::number).unwrap_or_else(Value::Error);
|
||||
let na = Value::Error(CellError::NA);
|
||||
for s in fit.se.iter().rev() {
|
||||
data.push(Value::number(*s));
|
||||
}
|
||||
data.push(fit.se_b.map(Value::number).unwrap_or(na.clone()));
|
||||
let rows: [[Value; 2]; 3] =
|
||||
[[rv(fit.r2), rv(fit.sey)], [rv(fit.f), Value::number(fit.df)], [Value::number(fit.ssreg), Value::number(fit.ssresid)]];
|
||||
for [p, q] in rows {
|
||||
data.push(p);
|
||||
if cols > 1 {
|
||||
data.push(q);
|
||||
}
|
||||
for _ in 2..cols {
|
||||
data.push(na.clone());
|
||||
}
|
||||
}
|
||||
return array_val(5, cols, data);
|
||||
}
|
||||
array_val(1, cols, data)
|
||||
}
|
||||
|
||||
pub(super) fn trend_impl(a: &[Arg], log: bool) -> R<Value> {
|
||||
let mut xy = parse_xy(&arg(a, 0)?.value, opt_value(a, 1))?;
|
||||
if log {
|
||||
if xy.y.iter().any(|&v| v <= 0.0) {
|
||||
return Err(CellError::Num);
|
||||
}
|
||||
xy.y.iter_mut().for_each(|v| *v = v.ln());
|
||||
}
|
||||
let constant = const_arg(a, 3)?;
|
||||
let fit = regress(&xy.y, &xy.x, xy.k, constant)?;
|
||||
let predict = |row: &[f64]| -> Value {
|
||||
let v = fit.b + fit.coef.iter().zip(row).map(|(c, x)| c * x).sum::<f64>();
|
||||
Value::number(if log { v.exp() } else { v })
|
||||
};
|
||||
match opt_value(a, 2) {
|
||||
None => {
|
||||
let data: Vec<Value> = xy.x.iter().map(|r| predict(r)).collect();
|
||||
array_val(xy.y_rows, xy.y_cols, data)
|
||||
}
|
||||
Some(nv) => {
|
||||
let na = as_array(nv);
|
||||
if xy.k == 1 {
|
||||
let data = na.data.iter().map(|v| strict_num(v).map(|x| predict(&[x]))).collect::<R<Vec<_>>>()?;
|
||||
array_val(na.rows, na.cols, data)
|
||||
} else if !xy.y_is_row && na.cols == xy.k {
|
||||
let mut data = Vec::with_capacity(na.rows);
|
||||
for r in 0..na.rows {
|
||||
let row = (0..na.cols).map(|c| strict_num(na.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?;
|
||||
data.push(predict(&row));
|
||||
}
|
||||
array_val(na.rows, 1, data)
|
||||
} else if xy.y_is_row && na.rows == xy.k {
|
||||
let mut data = Vec::with_capacity(na.cols);
|
||||
for c in 0..na.cols {
|
||||
let row = (0..na.rows).map(|r| strict_num(na.get(r, c).unwrap_or(&Value::Empty))).collect::<R<Vec<_>>>()?;
|
||||
data.push(predict(&row));
|
||||
}
|
||||
array_val(1, na.cols, data)
|
||||
} else {
|
||||
Err(CellError::Ref)
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user