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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
/// Answers "do these neighbors agree?" — a companion to a weighted mean.
///
/// A point estimate alone can't distinguish a neighborhood of tightly
/// agreeing observations from one where a few outliers and a pile of
/// unrelated noise happen to average out to the same number. `Dispersion`
/// reports the spread of neighbor outputs using the SAME weights as the
/// mean it accompanies (see [`crate::Neighbors::dispersion`] and
/// [`crate::Neighbors::gaussian_dispersion`]), so the two can never
/// describe different populations.
#[derive(Debug, Clone)]
pub struct Dispersion {
/// The weighted mean these statistics are computed about. Equal to
/// `weighted_mean()` (or `gaussian_weighted_mean(bandwidth)`) for the
/// same neighbor set.
pub mean: f64,
/// Weighted population variance of neighbor outputs about `mean`:
/// `Σw(x - mean)² / Σw`, using the same weights as `mean`. Zero when
/// every contributing neighbor agrees exactly.
pub variance: f64,
/// `sqrt(variance)`.
pub std_dev: f64,
/// Kish's effective sample size: `(Σw)² / Σw²`. This is
/// scale-invariant — k neighbors of any uniform weight report
/// `effective_n ≈ k` regardless of what that weight is, so k neighbors
/// at a uniform distance report ≈k no matter how far away that distance
/// is. It answers "how many roughly-independent observations back this
/// estimate", NOT "is there evidence near this query" — a large
/// `effective_n` built entirely from distant neighbors looks identical
/// to one built from close ones. Use `weight_sum` for a quantity that,
/// for a bounded kernel, distinguishes those.
pub effective_n: f64,
/// `Σw`: the total weight mass behind the estimate, using the same
/// weights as `mean`. For a kernel bounded in `[0, 1]` that decays with
/// distance (e.g. the Gaussian kernel behind
/// [`crate::Neighbors::gaussian_dispersion`]), this is a meaningful
/// "is there evidence near this query" signal: it shrinks toward 0 as
/// the query moves away from all training data. For the unbounded
/// inverse-distance kernel behind [`crate::Neighbors::dispersion`]
/// (`weight = instance_weight / distance`, which grows without bound as
/// distance approaches 0 and never reaches 0 as distance grows), it
/// does not have this property — prefer `effective_n` there.
pub weight_sum: f64,
}
impl Dispersion {
/// Compute dispersion statistics from `(weight, output)` pairs.
/// `pairs` must be non-empty — callers filter empty/degenerate cases
/// (no neighbors, or no neighbor contributing non-negligible weight)
/// before reaching here and return `None` instead.
pub(crate) fn from_weighted_pairs(pairs: &[(f64, f64)]) -> Self {
debug_assert!(
!pairs.is_empty(),
"Dispersion::from_weighted_pairs requires at least one pair"
);
let weight_sum: f64 = pairs.iter().map(|(w, _)| w).sum();
if weight_sum <= 0.0 {
// Degenerate: no positive weight to average over (e.g. every
// instance weight is zero or negative, which violates the
// crate's own "weight must be positive" contract but is
// enforced only by debug_assert in add_weighted). Fall back to
// the first pair alone, mirroring weighted_mean()'s own
// fallback in the equivalent situation.
return Dispersion {
mean: pairs[0].1,
variance: 0.0,
std_dev: 0.0,
effective_n: 1.0,
weight_sum,
};
}
let mean = pairs.iter().map(|(w, o)| w * o).sum::<f64>() / weight_sum;
let variance = pairs
.iter()
.map(|(w, o)| w * (o - mean) * (o - mean))
.sum::<f64>()
/ weight_sum;
// Kish's ESS, (Σw)²/Σw², is invariant to scaling every weight by the
// same positive constant — so compute it on weights normalized by
// their max instead of the raw weights. Mathematically identical,
// but avoids `w * w` overflowing to infinity (and the ratio
// collapsing to NaN) for legal large instance weights, since w² can
// overflow f64 long before w or Σw does.
let max_w = pairs.iter().fold(0.0_f64, |acc, (w, _)| acc.max(*w));
let effective_n = if max_w > 0.0 {
let scaled_sum: f64 = pairs.iter().map(|(w, _)| w / max_w).sum();
let scaled_sum_w2: f64 = pairs.iter().map(|(w, _)| (w / max_w).powi(2)).sum();
if scaled_sum_w2 > 0.0 {
scaled_sum * scaled_sum / scaled_sum_w2
} else {
0.0
}
} else {
0.0
};
Dispersion {
mean,
variance,
std_dev: variance.sqrt(),
effective_n,
weight_sum,
}
}
/// Standard error of `mean`: `sqrt(variance / effective_n)`.
///
/// Uses `effective_n` (Kish's effective sample size), never `weight_sum`
/// — this answers "how precisely is `mean` pinned down by roughly this
/// many independent-ish observations", which is `effective_n`'s job by
/// definition (see its doc comment). `weight_sum` is raw kernel mass;
/// dividing by it would conflate "how much evidence is nearby" with "how
/// many independent samples does that evidence represent", which for an
/// unbounded kernel (`dispersion()`) or a large uniform weight are not
/// the same number at all.
///
/// Caveat: this is a plug-in estimator computed entirely from the
/// neighbors it was built from, so it inherits their limits. With a
/// single neighbor (or several neighbors that all agree exactly),
/// `variance` is exactly zero by construction — there is nothing in the
/// neighbor set itself to measure spread against — so `standard_error()`
/// reports 0 (maximal confidence) rather than reflecting the real
/// uncertainty of estimating a mean from few observations. Callers that
/// need a noise floor for small neighborhoods should combine this with a
/// separate estimate (e.g. [`crate::Renegade::local_signal_variance`]'s
/// global comparison) rather than trusting `standard_error()` alone at
/// low `effective_n`.
pub fn standard_error(&self) -> f64 {
(self.variance / self.effective_n).sqrt()
}
/// Empirical-Bayes (James-Stein form) shrinkage of `mean` toward a
/// caller-supplied `prior`:
///
/// ```text
/// estimate = prior + λ · (mean − prior)
/// λ = signal_variance / (signal_variance + standard_error()²)
/// ```
///
/// `prior` is domain knowledge the crate has no way to know — a global
/// mean, a baseline rate, the model's own global prediction, whatever
/// the caller would fall back to with zero local evidence. `signal_variance`
/// is the between-neighborhood variance of the TRUE target: how much
/// legitimate local signal is there to trust, as opposed to noise.
/// [`crate::Renegade::local_signal_variance`] estimates it per-query from
/// the model's own data; see its docs for why a single global constant
/// is the wrong shape for this — it means using a domain-specific prior.
///
/// λ is the fraction of the gap between `prior` and `mean` that survives:
/// λ → 1 as the local estimate gets more precise (`standard_error` → 0)
/// or the neighborhood carries more real signal (`signal_variance` grows)
/// — trust the local mean fully. λ → 0 as the local estimate gets noisier
/// or the neighborhood carries no more signal than chance would produce
/// — fall back to the prior.
///
/// `signal_variance` must be non-negative; NaN propagates, a negative
/// non-NaN value is clamped to 0 (fully shrink to the prior — treated as
/// "no local signal detected" rather than an error). When both
/// `signal_variance` and `standard_error()` are exactly zero — no signal
/// estimate AND no measured spread (e.g. a single exact-match neighbor)
/// — there is nothing to distinguish trusting the local mean from
/// trusting the prior; this degenerate case defaults to λ = 1 (trust the
/// local observation), matching how `Dispersion` itself treats a single
/// pair as its own whole population. An infinite `signal_variance` (with
/// a finite noise term) takes the λ → 1 limit exactly instead of the
/// `∞/∞` a direct division would produce; symmetrically an infinite
/// noise term (`standard_error() = ∞`) with finite `signal_variance`
/// takes the λ → 0 limit. `prior`, `mean`, and `standard_error` are NOT
/// specially handled beyond that: a non-finite `prior` or `mean` flows
/// into `estimate = prior + λ·(mean − prior)` via ordinary IEEE-754
/// arithmetic (e.g. an infinite `prior` typically yields an infinite or
/// NaN `estimate`, depending on `λ` and `mean`), it is not forced to NaN.
pub fn shrink_toward(&self, prior: f64, signal_variance: f64) -> Shrinkage {
let signal_variance = if signal_variance.is_nan() {
signal_variance
} else {
signal_variance.max(0.0)
};
let standard_error = self.standard_error();
let noise_variance = standard_error * standard_error;
let lambda = if signal_variance.is_nan() || noise_variance.is_nan() {
f64::NAN
} else if signal_variance.is_infinite() && noise_variance.is_infinite() {
// ∞/∞: genuinely indeterminate, no limit to take.
f64::NAN
} else if signal_variance.is_infinite() {
1.0
} else if noise_variance.is_infinite() {
0.0
} else {
// Both finite. Normalize by the larger of the two before
// summing, rather than computing `signal_variance / (signal_variance
// + noise_variance)` directly: two individually-representable
// values (e.g. both near f64::MAX) can sum to +Infinity, which
// would silently zero out a ratio that should land near 0.5.
// Dividing both terms by their max first keeps the sum <= 2.0.
let scale = signal_variance.max(noise_variance);
if scale > 0.0 {
let sv = signal_variance / scale;
let nv = noise_variance / scale;
(sv / (sv + nv)).clamp(0.0, 1.0)
} else {
// Both exactly 0.
1.0
}
};
Shrinkage {
estimate: prior + lambda * (self.mean - prior),
lambda,
standard_error,
}
}
}
/// Result of [`Dispersion::shrink_toward`]: a local mean pulled toward a
/// prior by an amount that depends on how much the local evidence is worth
/// trusting.
#[derive(Debug, Clone)]
pub struct Shrinkage {
/// `prior + lambda * (mean - prior)`.
pub estimate: f64,
/// Shrinkage weight in `[0, 1]`. `1.0` = fully trust the local mean,
/// `0.0` = fully fall back to the prior.
pub lambda: f64,
/// `Dispersion::standard_error()` of the local mean being shrunk —
/// carried through so callers can see the precision behind `lambda`
/// without recomputing it.
pub standard_error: f64,
}