use bicmath_core::context::ExecContext;
use bicmath_core::contract::{
Args, CostClass, Example, FunctionDescriptor, Outcome, ParamDescriptor, SimpleFunction,
require_mode,
};
use bicmath_core::envelope::Exactness;
use bicmath_core::error::{EngineError, ErrorCode};
use bicmath_core::number::{Number, NumericMode};
use bicmath_core::schema::ValueSchema;
use bicmath_core::value::Value;
use num_bigint::BigInt;
use num_rational::BigRational;
use num_traits::Zero;
use std::sync::Arc;
use crate::common::*;
use crate::mathfn::sqrt;
fn moment_value(moment: &Moment) -> Result<Value, EngineError> {
match moment {
Moment::Exact(value) => Ok(rational_value(value.clone())),
Moment::Float(value) => float_value(*value),
}
}
fn count_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.count",
"statistics",
"1.0.0",
"Count",
"Count the observations in an array.",
)
.with_description(
"Returns the number of entries. Every entry must be a finite number; non-numeric or \
non-finite entries are rejected rather than silently skipped.",
)
.with_parameters(vec![ParamDescriptor::required(
"values",
"Observations to count.",
array_schema(any_number_schema()),
)])
.with_output(integer_schema(), "Number of observations.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#count")
.with_examples(vec![
Example::new(
"count three values",
example_args(&[("values", serde_json::json!([1, 2, 3]))]),
)
.with_value(parse_value(
serde_json::json!({"kind": "integer", "value": "3"}),
)),
])
}
fn invoke_count(args: &Args, _ctx: &ExecContext) -> Result<Outcome, EngineError> {
let values = args.array("values")?;
for (index, item) in values.iter().enumerate() {
match item {
Value::Number(number) => {
ensure_finite(number).map_err(|e| e.with_path(format!("values[{index}]")))?;
}
other => {
return Err(EngineError::malformed(format!(
"expected a number at values[{index}], found {}",
other.kind_name()
))
.with_path(format!("values[{index}]")));
}
}
}
Ok(Outcome::exact(Value::integer(BigInt::from(values.len()))))
}
fn mean_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.mean",
"statistics",
"1.0.0",
"Arithmetic mean",
"Exact arithmetic mean of a non-empty array.",
)
.with_description(
"Integer, rational, and decimal inputs produce an exact integer or rational mean. \
Float64 inputs require auto or scientific mode, are computed in binary64, and are \
marked approximate. An empty array is rejected with insufficient_observations.",
)
.with_parameters(vec![ParamDescriptor::required(
"values",
"Observations.",
array_schema(any_number_schema()),
)])
.with_output(any_number_schema(), "Arithmetic mean.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#mean")
.with_examples(vec![
Example::new(
"exact integer mean",
example_args(&[("values", serde_json::json!([1, 2, 3]))]),
)
.with_value(parse_value(
serde_json::json!({"kind": "integer", "value": "2"}),
)),
Example::new(
"exact rational mean",
example_args(&[("values", serde_json::json!([1, 2]))]),
)
.with_value(parse_value(
serde_json::json!({"kind": "rational", "numerator": "3", "denominator": "2"}),
)),
Example::new(
"empty array",
example_args(&[("values", serde_json::json!([]))]),
)
.with_error(ErrorCode::InsufficientObservations),
])
}
fn invoke_mean(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let series = classify_series(args, "values", ctx)?;
let moment = mean_moment(&series)?;
match moment {
Moment::Exact(value) => Ok(Outcome::exact(rational_value(value))),
Moment::Float(value) => Ok(Outcome::approximate(float_value(value)?)),
}
}
fn weighted_mean_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.weighted_mean",
"statistics",
"1.0.0",
"Weighted mean",
"Weighted arithmetic mean with frequency or reliability weights.",
)
.with_description(
"Frequency weights must be non-negative integers and count repeated observations; \
reliability weights must be non-negative numbers. Weights may be zero but the total \
weight must be positive. Negative weights, length mismatches, and a zero total weight \
are rejected. Exact inputs produce an exact rational result; float64 inputs require \
auto or scientific mode and are marked approximate.",
)
.with_parameters(vec![
ParamDescriptor::required("values", "Observations.", array_schema(any_number_schema())),
ParamDescriptor::required(
"weights",
"Non-negative weights, one per observation.",
array_schema(any_number_schema()),
),
ParamDescriptor::optional(
"weight_type",
"Weight semantics: frequency (default) or reliability.",
bicmath_core::schema::ValueSchema::Enum {
variants: vec!["frequency".to_string(), "reliability".to_string()],
},
),
])
.with_output(any_number_schema(), "Weighted mean.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#weighted_mean")
.with_examples(vec![
Example::new(
"frequency weights",
example_args(&[
("values", serde_json::json!([1, 2])),
("weights", serde_json::json!([1, 1])),
]),
)
.with_value(parse_value(
serde_json::json!({"kind": "rational", "numerator": "3", "denominator": "2"}),
)),
Example::new(
"negative weight",
example_args(&[
("values", serde_json::json!([1, 2])),
("weights", serde_json::json!([1, -1])),
]),
)
.with_error(ErrorCode::DomainViolation),
])
}
fn invoke_weighted_mean(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let weight_type = parse_choice(
args,
"weight_type",
"frequency",
&["frequency", "reliability"],
)?;
let values = collect_numbers(args, "values")?;
let weights = collect_numbers(args, "weights")?;
if values.len() != weights.len() {
return Err(EngineError::malformed(format!(
"values and weights must have the same length ({} vs {})",
values.len(),
weights.len()
)));
}
if values.is_empty() {
return Err(insufficient(
"weighted_mean requires at least one observation",
));
}
for (index, weight) in weights.iter().enumerate() {
if weight.is_negative() {
return Err(EngineError::domain("weights must be non-negative")
.with_path(format!("weights[{index}]")));
}
if weight_type == "frequency" {
let integral = match weight {
Number::Integer(_) => true,
Number::Rational(value) => value.is_integer(),
Number::Decimal(value) => value.to_bigint_if_integral().is_some(),
Number::Float64(value) => value.get().fract() == 0.0,
};
if !integral {
return Err(
EngineError::domain("frequency weights must be non-negative integers")
.with_path(format!("weights[{index}]")),
);
}
}
}
let has_float = values.iter().any(Number::is_float) || weights.iter().any(Number::is_float);
if has_float {
if ctx.numeric.mode == NumericMode::Exact {
return Err(EngineError::new(
ErrorCode::UnsupportedNumericMode,
"weighted_mean with float64 values or weights requires auto or scientific mode",
));
}
let mut numerator = 0.0f64;
let mut numerator_compensation = 0.0f64;
let mut denominator = 0.0f64;
let mut denominator_compensation = 0.0f64;
for (value, weight) in values.iter().zip(weights.iter()) {
let value = number_to_f64(value)?;
let weight = number_to_f64(weight)?;
let term = value * weight;
let adjusted = term - numerator_compensation;
let next = numerator + adjusted;
numerator_compensation = (next - numerator) - adjusted;
numerator = next;
let adjusted = weight - denominator_compensation;
let next = denominator + adjusted;
denominator_compensation = (next - denominator) - adjusted;
denominator = next;
}
if denominator == 0.0 {
return Err(EngineError::domain(
"total weight must be strictly positive",
));
}
return Ok(Outcome::approximate(float_value(numerator / denominator)?));
}
let mut numerator = BigRational::zero();
let mut denominator = BigRational::zero();
for (value, weight) in values.iter().zip(weights.iter()) {
let value = value
.as_exact_rational()
.ok_or_else(|| EngineError::internal("exact value could not be converted"))?;
let weight = weight
.as_exact_rational()
.ok_or_else(|| EngineError::internal("exact weight could not be converted"))?;
numerator += &value * &weight;
denominator += weight;
}
if denominator.is_zero() {
return Err(EngineError::domain(
"total weight must be strictly positive",
));
}
Ok(Outcome::exact(rational_value(numerator / denominator)))
}
fn median_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.median",
"statistics",
"1.0.0",
"Median",
"Median of a non-empty array.",
)
.with_description(
"The values are sorted numerically without modifying the input. An odd count returns \
the middle value; an even count returns the exact mean of the two middle values. \
Exact inputs stay exact (median([1, 2, 10, 100]) = 6); float64 inputs require auto \
or scientific mode.",
)
.with_parameters(vec![ParamDescriptor::required(
"values",
"Observations.",
array_schema(any_number_schema()),
)])
.with_output(any_number_schema(), "Median value.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#median")
.with_examples(vec![
Example::new(
"even count",
example_args(&[("values", serde_json::json!([1, 2, 10, 100]))]),
)
.with_value(parse_value(
serde_json::json!({"kind": "integer", "value": "6"}),
)),
Example::new(
"odd count",
example_args(&[("values", serde_json::json!([2, 10, 30]))]),
)
.with_value(parse_value(
serde_json::json!({"kind": "integer", "value": "10"}),
)),
])
}
fn invoke_median(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let series = classify_series(args, "values", ctx)?;
let (value, exactness) = median_number(&series)?;
Ok(Outcome::new(Value::Number(value), exactness))
}
fn mode_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.mode",
"statistics",
"1.0.0",
"Mode",
"All tied modes with their counts.",
)
.with_description(
"Returns every value that occurs most often, sorted ascending, with the matching \
counts. Ties are reported in full; there is no arbitrary winner. The method label is \
always all_modes.",
)
.with_parameters(vec![ParamDescriptor::required(
"values",
"Observations.",
array_schema(any_number_schema()),
)])
.with_output(
record_schema(
vec![
field("modes", array_schema(any_number_schema())),
field("counts", array_schema(integer_schema())),
field("method", text_schema()),
],
false,
),
"Record with modes, counts, and the method label all_modes.",
)
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#mode")
.with_examples(vec![
Example::new(
"single mode",
example_args(&[("values", serde_json::json!([1, 2, 2, 3]))]),
)
.with_value(parse_value(serde_json::json!({
"modes": [2],
"counts": [2],
"method": "all_modes"
}))),
])
}
fn invoke_mode(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let series = classify_series(args, "values", ctx)?;
let (modes, counts) = mode_numbers(&series)?;
let value = record(vec![
(
"modes",
array_value(modes.into_iter().map(number_value).collect()),
),
(
"counts",
array_value(
counts
.into_iter()
.map(|count| Value::integer(BigInt::from(count)))
.collect(),
),
),
("method", text("all_modes")),
]);
let exactness = if series.is_float() {
Exactness::Approximate
} else {
Exactness::Exact
};
Ok(Outcome::new(value, exactness))
}
fn minmax_descriptor(id: &str, title: &str, summary: &str, minimum: bool) -> FunctionDescriptor {
FunctionDescriptor::new(id, "statistics", "1.0.0", title, summary)
.with_description(
"Returns the smallest/largest observation without modifying the input. An empty \
array is rejected with insufficient_observations. Exact inputs return the exact \
selected value; float64 inputs require auto or scientific mode.",
)
.with_parameters(vec![ParamDescriptor::required(
"values",
"Observations.",
array_schema(any_number_schema()),
)])
.with_output(any_number_schema(), "Selected extreme value.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#min-max")
.with_examples(vec![
Example::new(
"extreme value",
example_args(&[("values", serde_json::json!([3, 1, 2]))]),
)
.with_value(parse_value(serde_json::json!({
"kind": "integer",
"value": if minimum { "1" } else { "3" }
}))),
])
}
fn invoke_min(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
select_minmax(args, ctx, true)
}
fn invoke_max(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
select_minmax(args, ctx, false)
}
fn select_minmax(args: &Args, ctx: &ExecContext, minimum: bool) -> Result<Outcome, EngineError> {
let series = classify_series(args, "values", ctx)?;
let (value, exactness) = extreme_number(&series, minimum)?;
Ok(Outcome::new(Value::Number(value), exactness))
}
fn quantile_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.quantile",
"statistics",
"1.0.0",
"Quantile",
"Quantile with an explicit interpolation convention.",
)
.with_description(
"Let h = (n - 1) * q. linear (default, R-7) returns \
x[floor(h)] + (h - floor(h)) * (x[floor(h)+1] - x[floor(h)]); lower returns \
x[floor(h)]; higher returns x[ceil(h)]; midpoint returns the mean of those two order \
statistics; nearest returns x[round_half_to_even(h)]. q must be in [0, 1]. Exact \
inputs stay exact; float64 inputs require auto or scientific mode.",
)
.with_parameters(vec![
ParamDescriptor::required("values", "Observations.", array_schema(any_number_schema())),
ParamDescriptor::required("q", "Quantile level in [0, 1].", any_number_schema()),
ParamDescriptor::optional(
"method",
"Interpolation method: linear (default), lower, higher, midpoint, nearest.",
bicmath_core::schema::ValueSchema::Enum {
variants: vec![
"linear".to_string(),
"lower".to_string(),
"higher".to_string(),
"midpoint".to_string(),
"nearest".to_string(),
],
},
),
])
.with_output(any_number_schema(), "Quantile value.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#quantile")
.with_examples(vec![
Example::new(
"linear interpolation",
example_args(&[
("values", serde_json::json!([1, 2, 3, 4])),
("q", serde_json::json!(0.5)),
]),
)
.with_value(parse_value(
serde_json::json!({"kind": "rational", "numerator": "5", "denominator": "2"}),
)),
Example::new(
"out of range q",
example_args(&[
("values", serde_json::json!([1, 2, 3])),
("q", serde_json::json!(1.5)),
]),
)
.with_error(ErrorCode::DomainViolation),
])
}
fn invoke_quantile(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let method = parse_choice(
args,
"method",
"linear",
&["linear", "lower", "higher", "midpoint", "nearest"],
)?;
let q = args.number("q")?.clone();
let series = classify_series(args, "values", ctx)?;
let (value, exactness) = quantile_number(&series, &q, &method)?;
Ok(Outcome::new(Value::Number(value), exactness))
}
fn variance_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.variance",
"statistics",
"1.0.0",
"Variance",
"Sample or population variance with explicit delta degrees of freedom.",
)
.with_description(
"Computes sum((x - mean)^2) / (n - ddof). ddof must be a non-negative integer and n \
must be strictly greater than ddof, otherwise insufficient_observations is returned. \
Exact inputs use exact rational arithmetic (variance([1, 2, 3], ddof=0) = 2/3); \
float64 inputs use Welford's stable one-pass algorithm and require auto or scientific \
mode.",
)
.with_parameters(vec![
ParamDescriptor::required("values", "Observations.", array_schema(any_number_schema())),
ParamDescriptor::required(
"ddof",
"Delta degrees of freedom: 0 for the population variance, 1 for the sample variance.",
integer_schema(),
),
])
.with_output(any_number_schema(), "Variance.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#variance")
.with_examples(vec![
Example::new(
"population variance",
example_args(&[
("values", serde_json::json!([1, 2, 3])),
("ddof", serde_json::json!(0)),
]),
)
.with_value(parse_value(
serde_json::json!({"kind": "rational", "numerator": "2", "denominator": "3"}),
)),
Example::new(
"sample variance",
example_args(&[
("values", serde_json::json!([1, 2, 3])),
("ddof", serde_json::json!(1)),
]),
)
.with_value(parse_value(
serde_json::json!({"kind": "integer", "value": "1"}),
)),
Example::new(
"too few observations",
example_args(&[
("values", serde_json::json!([1])),
("ddof", serde_json::json!(1)),
]),
)
.with_error(ErrorCode::InsufficientObservations),
])
}
fn invoke_variance(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let ddof = non_negative_u64(&args.integer("ddof")?, "ddof")?;
let series = classify_series(args, "values", ctx)?;
let moment = variance_moment(&series, ddof)?;
match moment {
Moment::Exact(value) => Ok(Outcome::exact(rational_value(value))),
Moment::Float(value) => Ok(Outcome::approximate(float_value(value)?)),
}
}
fn stddev_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.stddev",
"statistics",
"1.0.0",
"Standard deviation",
"Square root of the variance with explicit exactness handling.",
)
.with_description(
"Computes sqrt(variance(values, ddof)). When the variance is a perfect square in the \
selected representation the result is exact. Otherwise exact mode returns \
unsupported_numeric_mode, auto mode returns a decimal approximation marked \
approximate, and scientific mode returns float64. Float64 inputs require auto or \
scientific mode.",
)
.with_parameters(vec![
ParamDescriptor::required("values", "Observations.", array_schema(any_number_schema())),
ParamDescriptor::required(
"ddof",
"Delta degrees of freedom: 0 for the population standard deviation, 1 for the sample.",
integer_schema(),
),
])
.with_output(any_number_schema(), "Standard deviation.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#stddev")
.with_examples(vec![
Example::new(
"exact perfect square",
example_args(&[
("values", serde_json::json!([1, 2, 3])),
("ddof", serde_json::json!(1)),
]),
)
.with_value(parse_value(
serde_json::json!({"kind": "integer", "value": "1"}),
)),
Example::new(
"non-square variance in exact mode",
example_args(&[
("values", serde_json::json!([1, 2, 3])),
("ddof", serde_json::json!(0)),
]),
)
.with_error(ErrorCode::UnsupportedNumericMode),
])
}
fn invoke_stddev(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let ddof = non_negative_u64(&args.integer("ddof")?, "ddof")?;
let series = classify_series(args, "values", ctx)?;
let moment = variance_moment(&series, ddof)?;
let (value, exactness) = stddev_moment(&moment, ctx.numeric.mode)?;
Ok(Outcome::new(Value::Number(value), exactness))
}
fn covariance_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.covariance",
"statistics",
"1.0.0",
"Covariance",
"Sample or population covariance of two paired series.",
)
.with_description(
"Computes sum((x - mean_x) * (y - mean_y)) / (n - ddof) for paired observations. The \
two arrays must have equal length; n must be strictly greater than ddof. Exact inputs \
use exact rational arithmetic; float64 inputs use a stable two-pass centered \
algorithm and require auto or scientific mode.",
)
.with_parameters(vec![
ParamDescriptor::required("xs", "First series.", array_schema(any_number_schema())),
ParamDescriptor::required(
"ys",
"Second series, paired with xs.",
array_schema(any_number_schema()),
),
ParamDescriptor::required("ddof", "Delta degrees of freedom.", integer_schema()),
])
.with_output(any_number_schema(), "Covariance.")
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#covariance")
.with_examples(vec![
Example::new(
"covariance of identical series",
example_args(&[
("xs", serde_json::json!([1, 2, 3])),
("ys", serde_json::json!([1, 2, 3])),
("ddof", serde_json::json!(0)),
]),
)
.with_value(parse_value(
serde_json::json!({"kind": "rational", "numerator": "2", "denominator": "3"}),
)),
Example::new(
"length mismatch",
example_args(&[
("xs", serde_json::json!([1, 2, 3])),
("ys", serde_json::json!([1, 2])),
("ddof", serde_json::json!(0)),
]),
)
.with_error(ErrorCode::MalformedInput),
])
}
fn invoke_covariance(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let ddof = non_negative_u64(&args.integer("ddof")?, "ddof")?;
let xs = classify_series(args, "xs", ctx)?;
let ys = classify_series(args, "ys", ctx)?;
if xs.len() != ys.len() {
return Err(EngineError::malformed(format!(
"xs and ys must have the same length ({} vs {})",
xs.len(),
ys.len()
)));
}
let moment = covariance_moment(&xs, &ys, ddof)?;
match moment {
Moment::Exact(value) => Ok(Outcome::exact(rational_value(value))),
Moment::Float(value) => Ok(Outcome::approximate(float_value(value)?)),
}
}
fn correlation_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.correlation",
"statistics",
"1.0.0",
"Pearson correlation",
"Pearson product-moment correlation coefficient.",
)
.with_description(
"Returns the Pearson correlation r = cov(x, y) / (sd(x) * sd(y)) as float64. If either \
series has zero variance the result is undefined and a domain_violation is returned \
instead of NaN. At least two paired observations are required. The output is always \
approximate, so exact mode is not supported.",
)
.with_parameters(vec![
ParamDescriptor::required("xs", "First series.", array_schema(any_number_schema())),
ParamDescriptor::required(
"ys",
"Second series, paired with xs.",
array_schema(any_number_schema()),
),
])
.with_output(float64_schema(), "Pearson correlation coefficient.")
.with_modes(inferential_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#correlation")
.with_examples(vec![
Example::new(
"perfect positive correlation",
example_args(&[
("xs", serde_json::json!([1, 2, 3])),
("ys", serde_json::json!([2, 4, 6])),
]),
)
.with_value(parse_value(
serde_json::json!({"kind": "float64", "value": "1"}),
)),
Example::new(
"zero variance",
example_args(&[
("xs", serde_json::json!([1, 1, 1])),
("ys", serde_json::json!([1, 2, 3])),
]),
)
.with_error(ErrorCode::DomainViolation),
])
}
fn invoke_correlation(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
require_mode(ctx, &inferential_modes(), "statistics.correlation")?;
let xs = collect_numbers(args, "xs")?;
let ys = collect_numbers(args, "ys")?;
if xs.len() != ys.len() {
return Err(EngineError::malformed(format!(
"xs and ys must have the same length ({} vs {})",
xs.len(),
ys.len()
)));
}
if xs.len() < 2 {
return Err(insufficient(
"correlation requires at least two paired observations",
));
}
let xs: Vec<f64> = xs
.iter()
.enumerate()
.map(|(index, value)| number_to_f64(value).map_err(|e| e.with_path(format!("xs[{index}]"))))
.collect::<Result<_, _>>()?;
let ys: Vec<f64> = ys
.iter()
.enumerate()
.map(|(index, value)| number_to_f64(value).map_err(|e| e.with_path(format!("ys[{index}]"))))
.collect::<Result<_, _>>()?;
let n = xs.len();
let mean_x = mean_f64(&xs);
let mean_y = mean_f64(&ys);
let mut sxx = 0.0f64;
let mut syy = 0.0f64;
let mut sxy = 0.0f64;
for index in 0..n {
let dx = xs[index] - mean_x;
let dy = ys[index] - mean_y;
sxx += dx * dx;
syy += dy * dy;
sxy += dx * dy;
}
if sxx <= 0.0 || syy <= 0.0 {
return Err(EngineError::domain(
"correlation is undefined when either series has zero variance",
));
}
let denominator = sqrt(sxx * syy);
if denominator <= 0.0 {
return Err(EngineError::domain(
"correlation is undefined when either series has zero variance",
));
}
let r = (sxy / denominator).clamp(-1.0, 1.0);
Ok(Outcome::approximate(float_value(r)?))
}
fn summary_descriptor() -> FunctionDescriptor {
FunctionDescriptor::new(
"statistics.summary",
"statistics",
"1.0.0",
"Summary statistics",
"One-call descriptive summary with quantiles.",
)
.with_description(
"Returns count, mean, min, max, median, sample variance (ddof = 1), sample standard \
deviation, and the linear quantiles 0.25, 0.5, and 0.75. At least two observations \
are required. Exact inputs stay exact wherever the representation allows; the \
standard deviation follows the stddev exactness rules.",
)
.with_parameters(vec![ParamDescriptor::required(
"values",
"Observations.",
array_schema(any_number_schema()),
)])
.with_output(
record_schema(
vec![
field("count", integer_schema()),
field("mean", any_number_schema()),
field("min", any_number_schema()),
field("max", any_number_schema()),
field("median", any_number_schema()),
field("variance", any_number_schema()),
field("stddev", any_number_schema()),
field("quantiles", ValueSchema::Any),
],
false,
),
"Record with count, mean, min, max, median, variance, stddev, and quantiles.",
)
.with_modes(all_modes())
.with_cost(CostClass::Linear)
.with_method_ref("docs/methods/statistics.md#summary")
.with_examples(vec![
Example::new(
"too few observations",
example_args(&[("values", serde_json::json!([1]))]),
)
.with_error(ErrorCode::InsufficientObservations),
])
}
fn invoke_summary(args: &Args, ctx: &ExecContext) -> Result<Outcome, EngineError> {
let series = classify_series(args, "values", ctx)?;
if series.len() < 2 {
return Err(insufficient(
"summary requires at least two observations because variance uses ddof = 1",
));
}
let mean = mean_moment(&series)?;
let variance = variance_moment(&series, 1)?;
let (stddev, stddev_exactness) = stddev_moment(&variance, ctx.numeric.mode)?;
let (minimum, _) = extreme_number(&series, true)?;
let (maximum, _) = extreme_number(&series, false)?;
let (median, _) = median_number(&series)?;
let quarter = rational_to_number(BigRational::new(BigInt::from(1), BigInt::from(4)));
let half = rational_to_number(BigRational::new(BigInt::from(1), BigInt::from(2)));
let three_quarters = rational_to_number(BigRational::new(BigInt::from(3), BigInt::from(4)));
let (q25, _) = quantile_number(&series, &quarter, "linear")?;
let (q50, _) = quantile_number(&series, &half, "linear")?;
let (q75, _) = quantile_number(&series, &three_quarters, "linear")?;
let value = record(vec![
("count", integer_value(series.len() as u64)),
("mean", moment_value(&mean)?),
("min", Value::Number(minimum)),
("max", Value::Number(maximum)),
("median", Value::Number(median)),
("variance", moment_value(&variance)?),
("stddev", Value::Number(stddev)),
(
"quantiles",
record(vec![
("0.25", Value::Number(q25)),
("0.5", Value::Number(q50)),
("0.75", Value::Number(q75)),
]),
),
]);
let exactness = Exactness::Exact.combine(stddev_exactness);
Ok(Outcome::new(value, exactness))
}
pub fn functions() -> Vec<Arc<dyn bicmath_core::contract::Function>> {
vec![
SimpleFunction::arc(count_descriptor(), invoke_count),
SimpleFunction::arc(mean_descriptor(), invoke_mean),
SimpleFunction::arc(weighted_mean_descriptor(), invoke_weighted_mean),
SimpleFunction::arc(median_descriptor(), invoke_median),
SimpleFunction::arc(mode_descriptor(), invoke_mode),
SimpleFunction::arc(
minmax_descriptor(
"statistics.min",
"Minimum",
"Smallest observation in a non-empty array.",
true,
),
invoke_min,
),
SimpleFunction::arc(
minmax_descriptor(
"statistics.max",
"Maximum",
"Largest observation in a non-empty array.",
false,
),
invoke_max,
),
SimpleFunction::arc(quantile_descriptor(), invoke_quantile),
SimpleFunction::arc(variance_descriptor(), invoke_variance),
SimpleFunction::arc(stddev_descriptor(), invoke_stddev),
SimpleFunction::arc(covariance_descriptor(), invoke_covariance),
SimpleFunction::arc(correlation_descriptor(), invoke_correlation),
SimpleFunction::arc(summary_descriptor(), invoke_summary),
]
}
#[cfg(test)]
mod tests {
use super::*;
use std::collections::BTreeMap;
fn call(id: &str, raw: serde_json::Value) -> Result<Outcome, EngineError> {
call_with_mode(id, raw, NumericMode::Auto)
}
fn call_with_mode(
id: &str,
raw: serde_json::Value,
mode: NumericMode,
) -> Result<Outcome, EngineError> {
let module = crate::module();
let function = module
.functions
.iter()
.find(|f| f.descriptor().id == id)
.expect("function exists");
let ctx = ExecContext {
numeric: bicmath_core::number::NumericContext {
mode,
..bicmath_core::number::NumericContext::default()
},
..ExecContext::conservative()
};
let args_json = raw.as_object().expect("object args");
let mut values = BTreeMap::new();
for (name, value) in args_json {
let param = function
.descriptor()
.parameter(name)
.expect("parameter exists");
values.insert(
name.clone(),
param
.schema
.coerce(value, name, &ctx.limits, true)
.expect("argument coerces"),
);
}
function.invoke(&Args::new(values), &ctx)
}
fn text_of(outcome: &Outcome) -> String {
match &outcome.value {
Value::Number(number) => number.to_string(),
other => panic!("expected number, got {other:?}"),
}
}
#[test]
fn exact_descriptive_values() {
assert_eq!(
text_of(&call("statistics.mean", serde_json::json!({"values": [1, 2, 3]})).unwrap()),
"2"
);
assert_eq!(
text_of(&call("statistics.mean", serde_json::json!({"values": [1, 2]})).unwrap()),
"3/2"
);
assert_eq!(
text_of(
&call(
"statistics.median",
serde_json::json!({"values": [2, 10, 30]})
)
.unwrap()
),
"10"
);
assert_eq!(
text_of(
&call(
"statistics.median",
serde_json::json!({"values": [1, 2, 10, 100]})
)
.unwrap()
),
"6"
);
assert_eq!(
text_of(
&call(
"statistics.variance",
serde_json::json!({"values": [1, 2, 3], "ddof": 0})
)
.unwrap()
),
"2/3"
);
assert_eq!(
text_of(
&call(
"statistics.variance",
serde_json::json!({"values": [1, 2, 3], "ddof": 1})
)
.unwrap()
),
"1"
);
}
#[test]
fn original_input_is_not_mutated() {
let raw = serde_json::json!({"values": [3, 1, 2]});
let outcome = call("statistics.median", raw.clone()).unwrap();
assert_eq!(text_of(&outcome), "2");
let args = raw.as_object().unwrap();
assert_eq!(args["values"], serde_json::json!([3, 1, 2]));
}
#[test]
fn float_inputs_are_approximate() {
let outcome = call(
"statistics.mean",
serde_json::json!({"values": [
{"kind": "float64", "value": "1"},
{"kind": "float64", "value": "2"}
]}),
)
.unwrap();
assert_eq!(outcome.exactness, Exactness::Approximate);
assert_eq!(text_of(&outcome), "1.5");
}
#[test]
fn exact_mode_rejects_float_inputs() {
let err = call_with_mode(
"statistics.mean",
serde_json::json!({"values": [{"kind": "float64", "value": "1"}]}),
NumericMode::Exact,
)
.unwrap_err();
assert_eq!(err.code, ErrorCode::UnsupportedNumericMode);
}
#[test]
fn non_numeric_entries_are_rejected() {
let module = crate::module();
let function = module
.functions
.iter()
.find(|f| f.descriptor().id == "statistics.mean")
.expect("function exists");
let mut values = BTreeMap::new();
values.insert(
"values".to_string(),
Value::Array(vec![
Value::Number(Number::Integer(BigInt::from(1))),
Value::text("not a number"),
]),
);
let err = function
.invoke(&Args::new(values), &ExecContext::conservative())
.unwrap_err();
assert_eq!(err.code, ErrorCode::MalformedInput);
}
#[test]
fn mode_returns_all_ties() {
let outcome = call(
"statistics.mode",
serde_json::json!({"values": [1, 2, 2, 3, 3]}),
)
.unwrap();
let Value::Record(fields) = &outcome.value else {
panic!("expected record");
};
assert_eq!(
fields.get("modes"),
Some(&Value::Array(vec![
Value::Number(Number::Integer(BigInt::from(2))),
Value::Number(Number::Integer(BigInt::from(3))),
]))
);
}
#[test]
fn weighted_mean_rejects_bad_weights() {
let err = call(
"statistics.weighted_mean",
serde_json::json!({"values": [1, 2], "weights": [1, -1]}),
)
.unwrap_err();
assert_eq!(err.code, ErrorCode::DomainViolation);
let err = call(
"statistics.weighted_mean",
serde_json::json!({"values": [1, 2], "weights": [0, 0]}),
)
.unwrap_err();
assert_eq!(err.code, ErrorCode::DomainViolation);
let err = call(
"statistics.weighted_mean",
serde_json::json!({"values": [1, 2], "weights": [1]}),
)
.unwrap_err();
assert_eq!(err.code, ErrorCode::MalformedInput);
}
#[test]
fn stddev_exact_and_approximate_paths() {
let exact = call(
"statistics.stddev",
serde_json::json!({"values": [1, 2, 3], "ddof": 1}),
)
.unwrap();
assert_eq!(exact.exactness, Exactness::Exact);
let err = call_with_mode(
"statistics.stddev",
serde_json::json!({"values": [1, 2, 3], "ddof": 0}),
NumericMode::Exact,
)
.unwrap_err();
assert_eq!(err.code, ErrorCode::UnsupportedNumericMode);
let approximate = call(
"statistics.stddev",
serde_json::json!({"values": [1, 2, 3], "ddof": 0}),
)
.unwrap();
assert_eq!(approximate.exactness, Exactness::Approximate);
}
#[test]
fn correlation_detects_zero_variance() {
let err = call(
"statistics.correlation",
serde_json::json!({"xs": [1, 1, 1], "ys": [1, 2, 3]}),
)
.unwrap_err();
assert_eq!(err.code, ErrorCode::DomainViolation);
}
#[test]
fn quantile_methods() {
let lower = call(
"statistics.quantile",
serde_json::json!({"values": [1, 2, 3, 4], "q": "0.5", "method": "lower"}),
)
.unwrap();
assert_eq!(text_of(&lower), "2");
let higher = call(
"statistics.quantile",
serde_json::json!({"values": [1, 2, 3, 4], "q": "0.5", "method": "higher"}),
)
.unwrap();
assert_eq!(text_of(&higher), "3");
let nearest = call(
"statistics.quantile",
serde_json::json!({"values": [1, 2, 3, 4], "q": "0.5", "method": "nearest"}),
)
.unwrap();
assert_eq!(text_of(&nearest), "3");
}
#[test]
fn summary_exact_record() {
let outcome = call(
"statistics.summary",
serde_json::json!({"values": [1, 2, 3]}),
)
.unwrap();
assert_eq!(outcome.exactness, Exactness::Exact);
let Value::Record(fields) = &outcome.value else {
panic!("expected record");
};
assert_eq!(
fields.get("variance"),
Some(&Value::Number(Number::Integer(BigInt::from(1))))
);
}
}