mirror of
https://github.com/storytold/gridcraft.git
synced 2026-10-31 14:17:12 -06:00
732 lines
21 KiB
Rust
732 lines
21 KiB
Rust
use super::descriptive::{arg_nums, arg_pairs, mean, var_s};
|
||
use super::{
|
||
Arg, C, CellError, R, Value, arg, as_array, boolean, erf, gamma, has, invert, ln_beta, ln_gamma, norm_cdf, norm_inv, norm_pdf, num, num_val,
|
||
reg_gamma_p, reg_gamma_q, reg_inc_beta,
|
||
};
|
||
|
||
// ---------------------------------------------------------------------------------------------
|
||
// 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 {
|
||
let z = v / 2.0;
|
||
let log_gamma_ratio = if z >= 16.0 {
|
||
let inv_z = 1.0 / z;
|
||
0.5 * z.ln() - inv_z / 8.0 + inv_z.powi(3) / 192.0 - inv_z.powi(5) / 640.0 + 17.0 * inv_z.powi(7) / 14_336.0
|
||
} else {
|
||
ln_gamma(z + 0.5) - ln_gamma(z)
|
||
};
|
||
(log_gamma_ratio - 0.5 * (v.ln() + std::f64::consts::PI.ln()) - (v + 1.0) / 2.0 * (t * t / v).ln_1p()).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;
|
||
}
|
||
if x == 1.0 && d1 == d2 {
|
||
return 0.5;
|
||
}
|
||
if d1 >= 1e8 && d2 / d1 <= 1e-8 {
|
||
return reg_gamma_q(d2 / 2.0, d2 / (2.0 * (x + d2 / d1)));
|
||
}
|
||
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;
|
||
}
|
||
if d1 >= 1e8 && d2 / d1 <= 1e-8 {
|
||
// With a much larger numerator degree, the scaled beta tail converges to a gamma tail.
|
||
return reg_gamma_p(d2 / 2.0, d2 / (2.0 * (x + d2 / d1)));
|
||
}
|
||
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()
|
||
}
|
||
|
||
const MAX_BINOM_TRIALS: f64 = 2_147_483_646.0;
|
||
|
||
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 = num(a, 1)?.trunc();
|
||
let cum = boolean(a, 2)?;
|
||
if v == 0.0 {
|
||
return Err(if cum { CellError::Num } else { CellError::Div0 });
|
||
}
|
||
if v < 0.0 {
|
||
return Err(CellError::Num);
|
||
}
|
||
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 {
|
||
return Err(CellError::Div0);
|
||
}
|
||
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);
|
||
}
|
||
let value = if p < 0.5 { invert(|x| f_cdf(x, d1, d2), p, 0.0, 1.0, 1e60) } else { invert(|x| f_upper(x, d1, d2), 1.0 - p, 0.0, 1.0, 1e60) };
|
||
num_val(value)
|
||
}
|
||
|
||
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 || n > MAX_BINOM_TRIALS || 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..=MAX_BINOM_TRIALS).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 { 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);
|
||
}
|
||
if x == 0.0 {
|
||
return num_val(0.0);
|
||
}
|
||
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))
|
||
}
|