1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
//! Streaming ordinary least squares: slope, intercept, and R² in one pass.
/// A least-squares line fit over a stream of `(x, y)` pairs: bivariate
/// Welford accumulation of the means, variances, and covariance, from which
/// the slope, intercept, coefficient of determination, and the standard error
/// of the mean response derive.
///
/// A mergeable monoid: partial results over chunks combine with
/// [`Fit::merge`] (Chan's parallel formula) into the same statistics a single
/// pass produces — host-side parallelism and streaming are compositions.
/// Pairs with a non-finite member are ignored, matching the gap convention.
#[derive(Debug, Clone, Copy)]
pub struct Fit {
count: u64,
mean_x: f64,
mean_y: f64,
/// Σ(x−x̄)² — the x sum of squares.
m2_x: f64,
/// Σ(y−ȳ)² — the y sum of squares.
m2_y: f64,
/// Σ(x−x̄)(y−ȳ) — the co-moment.
co: f64,
}
impl Fit {
/// An empty accumulator.
pub fn new() -> Fit {
Fit {
count: 0,
mean_x: 0.0,
mean_y: 0.0,
m2_x: 0.0,
m2_y: 0.0,
co: 0.0,
}
}
/// The fit of paired slices — [`Fit::add`] over `x.iter().zip(y)`.
///
/// # Panics
///
/// Panics if the slices have different lengths.
pub fn xy(x: &[f64], y: &[f64]) -> Fit {
assert_eq!(x.len(), y.len(), "Fit::xy requires slices of equal length");
let mut fit = Fit::new();
for (&x, &y) in x.iter().zip(y) {
fit.add(x, y);
}
fit
}
/// Accumulates one pair; pairs with a non-finite member are ignored.
pub fn add(&mut self, x: f64, y: f64) {
if !x.is_finite() || !y.is_finite() {
return;
}
self.count += 1;
let count = self.count as f64;
let delta_x = x - self.mean_x;
let delta_y = y - self.mean_y;
self.mean_x += delta_x / count;
self.mean_y += delta_y / count;
// The second factors use the updated means (Welford's trick), keeping
// every accumulator cancellation-free for offset data like timestamps.
self.m2_x += delta_x * (x - self.mean_x);
self.m2_y += delta_y * (y - self.mean_y);
self.co += delta_x * (y - self.mean_y);
}
/// Merges another accumulator into this one.
pub fn merge(&mut self, other: &Fit) {
if other.count == 0 {
return;
}
if self.count == 0 {
*self = *other;
return;
}
let total = self.count + other.count;
let weight = self.count as f64 * other.count as f64 / total as f64;
let delta_x = other.mean_x - self.mean_x;
let delta_y = other.mean_y - self.mean_y;
self.m2_x += other.m2_x + delta_x * delta_x * weight;
self.m2_y += other.m2_y + delta_y * delta_y * weight;
self.co += other.co + delta_x * delta_y * weight;
self.mean_x += delta_x * other.count as f64 / total as f64;
self.mean_y += delta_y * other.count as f64 / total as f64;
self.count = total;
}
/// The number of finite pairs seen.
pub fn count(&self) -> u64 {
self.count
}
/// The least-squares slope, or `None` while the fit is degenerate (fewer
/// than two pairs, or no x spread).
pub fn slope(&self) -> Option<f64> {
(self.count >= 2 && self.m2_x > 0.0).then(|| self.co / self.m2_x)
}
/// The least-squares intercept, or `None` while the fit is degenerate.
pub fn intercept(&self) -> Option<f64> {
self.slope().map(|slope| self.mean_y - slope * self.mean_x)
}
/// The fitted value at `x`, or `None` while the fit is degenerate.
pub fn predict(&self, x: f64) -> Option<f64> {
self.slope()
.map(|slope| self.mean_y + slope * (x - self.mean_x))
}
/// The coefficient of determination R², or `None` while the fit is
/// degenerate. Constant y fits perfectly: R² is 1.
pub fn r_squared(&self) -> Option<f64> {
self.slope()?;
if self.m2_y == 0.0 {
return Some(1.0);
}
Some((self.co * self.co / (self.m2_x * self.m2_y)).clamp(0.0, 1.0))
}
/// The standard error of the *mean response* at `x` — the half-width unit
/// of a confidence band around the fitted line (multiply by a quantile:
/// 1.96 approximates a 95% band for large samples). `None` while the fit
/// is degenerate or has no residual degrees of freedom (fewer than three
/// pairs).
pub fn standard_error(&self, x: f64) -> Option<f64> {
self.slope()?;
if self.count < 3 {
return None;
}
let residual = (self.m2_y - self.co * self.co / self.m2_x).max(0.0);
let variance = residual / (self.count - 2) as f64;
let offset = x - self.mean_x;
let leverage = 1.0 / self.count as f64 + offset * offset / self.m2_x;
Some((variance * leverage).sqrt())
}
}
impl Default for Fit {
/// The empty accumulator — identical to [`Fit::new`].
fn default() -> Fit {
Fit::new()
}
}
#[cfg(test)]
#[path = "tests/fit_tests.rs"]
mod tests;