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 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302
#![allow(missing_docs)]
use std::mem;
#[cfg(test)]
mod tests;
fn local_sort(v: &mut [f64]) {
v.sort_by(|x: &f64, y: &f64| x.total_cmp(y));
}
/// 一个 trait,它提供关于单变量数字样本集的简单描述性统计信息。
pub trait Stats {
/// 样本的总和。
///
/// Note: 这种方法牺牲了性能,准确性取决于 IEEE 754 算法保证。
/// 请参见以下网址的正确性证明:
/// ["Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates"][paper]
///
/// [paper]: https://www.cs.cmu.edu/~quake-papers/robust-arithmetic.ps
///
fn sum(&self) -> f64;
/// 样品的最小值。
fn min(&self) -> f64;
/// 样品的最大值。
fn max(&self) -> f64;
/// 样本的算术平均值 (average): 总和除以样本数。
///
/// 请参见:<https://en.wikipedia.org/wiki/Arithmetic_mean>
fn mean(&self) -> f64;
/// 样本中位数:将样本下半部分与样本上半部分分开的值。
/// 等于 `self.percentile(50.0)`。
///
/// 请参见:<https://en.wikipedia.org/wiki/Median>
fn median(&self) -> f64;
/// 样本方差:每个样本与样本均值之差的平方的偏差校正后的平均值。
/// 请注意,这是计算 _sample variance_ 而不是总体方差 (假定该方差是未知的)。
/// 因此,通过除以 `(n-1)` 而不是 `n`,它可以纠正如果我们计算总体方差时出现的 `(n-1)/n` 偏差。
///
///
/// 请参见:<https://en.wikipedia.org/wiki/Variance>
///
fn var(&self) -> f64;
/// 标准偏差:样本方差的平方根。
///
/// Note: 对于非正态分布,这不是一个可靠的统计数据。
/// 对于未知的发行版,最好使用 `median_abs_dev`。
///
/// 请参见:<https://en.wikipedia.org/wiki/Standard_deviation>
fn std_dev(&self) -> f64;
/// 标准偏差,以平均值的百分比表示。请参见 `std_dev` 和 `mean`。
///
/// Note: 对于非正态分布,这不是一个可靠的统计数据。
/// 对于未知的发行版,最好使用 `median_abs_dev_pct`。
fn std_dev_pct(&self) -> f64;
/// 每个样本与样本中位数的绝对偏差的标度中位数。
/// 这是样本可变性的可靠 (distribution-agnostic) 估计器。
/// 如果您不能假定样本是正态分布的,则优先使用 `std_dev`。
/// 请注意,常量 `1.4826` 对它进行了缩放,以使其可用作标准偏差的一致估计量。
///
///
/// 请参见:<https://en.wikipedia.org/wiki/Median_absolute_deviation>
fn median_abs_dev(&self) -> f64;
/// 中位数绝对偏差占中位数的百分比。请参见 `median_abs_dev` 和 `median`。
fn median_abs_dev_pct(&self) -> f64;
/// 百分位数: `self` 中的值所占的 `pct` 百分比以下的值。
/// 例如,percentile(95.0) 将返回值 `v`,以使 `self` 中 95% 的样本 `s` 满足 `s <= v`。
///
///
/// 通过最接近等级之间的线性插值计算。
///
/// 请参见:<https://en.wikipedia.org/wiki/Percentile>
fn percentile(&self, pct: f64) -> f64;
/// 样本四分位数:将样本分为四个相等组的三个值,每个组的数据为 1/4。
/// 中间值为中间值。
/// 请参见 `median` 和 `percentile`。
/// 与对 `percentile` 的 3 个调用相比,此函数可以更有效地计算 3 个四分位数,但在其他方面等效。
///
/// 另请参见:<https://en.wikipedia.org/wiki/Quartile>
fn quartiles(&self) -> (f64, f64, f64);
/// 四分位间距:第 25 个百分位数 (第 1 个四分位数) 和第 75 个百分位数 (第 3 个四分位数) 之间的差。
/// 请参见 `quartiles`。
///
/// 另请参见:<https://en.wikipedia.org/wiki/Interquartile_range>
fn iqr(&self) -> f64;
}
/// 样本集的所有汇总统计信息的提取集合。
#[derive(Debug, Clone, PartialEq, Copy)]
#[allow(missing_docs)]
pub struct Summary {
pub sum: f64,
pub min: f64,
pub max: f64,
pub mean: f64,
pub median: f64,
pub var: f64,
pub std_dev: f64,
pub std_dev_pct: f64,
pub median_abs_dev: f64,
pub median_abs_dev_pct: f64,
pub quartiles: (f64, f64, f64),
pub iqr: f64,
}
impl Summary {
/// 创建一个新的样本集总结。
pub fn new(samples: &[f64]) -> Summary {
Summary {
sum: samples.sum(),
min: samples.min(),
max: samples.max(),
mean: samples.mean(),
median: samples.median(),
var: samples.var(),
std_dev: samples.std_dev(),
std_dev_pct: samples.std_dev_pct(),
median_abs_dev: samples.median_abs_dev(),
median_abs_dev_pct: samples.median_abs_dev_pct(),
quartiles: samples.quartiles(),
iqr: samples.iqr(),
}
}
}
impl Stats for [f64] {
// FIXME #11059 处理 NaN,inf 和溢出
fn sum(&self) -> f64 {
let mut partials = vec![];
for &x in self {
let mut x = x;
let mut j = 0;
// 此内部循环将 `hi`/`lo` 总和应用于每个部分,以使部分总和的列表保持准确。
//
for i in 0..partials.len() {
let mut y: f64 = partials[i];
if x.abs() < y.abs() {
mem::swap(&mut x, &mut y);
}
// 四舍五入的 `x+y` 存储在 `hi` 中,四舍五入存储在 `lo` 中。
// `hi+lo` 等于 `x+y`。
let hi = x + y;
let lo = y - (hi - x);
if lo != 0.0 {
partials[j] = lo;
j += 1;
}
x = hi;
}
if j >= partials.len() {
partials.push(x);
} else {
partials[j] = x;
partials.truncate(j + 1);
}
}
let zero: f64 = 0.0;
partials.iter().fold(zero, |p, q| p + *q)
}
fn min(&self) -> f64 {
assert!(!self.is_empty());
self.iter().fold(self[0], |p, q| p.min(*q))
}
fn max(&self) -> f64 {
assert!(!self.is_empty());
self.iter().fold(self[0], |p, q| p.max(*q))
}
fn mean(&self) -> f64 {
assert!(!self.is_empty());
self.sum() / (self.len() as f64)
}
fn median(&self) -> f64 {
self.percentile(50_f64)
}
fn var(&self) -> f64 {
if self.len() < 2 {
0.0
} else {
let mean = self.mean();
let mut v: f64 = 0.0;
for s in self {
let x = *s - mean;
v += x * x;
}
// 注意,这 _应该是 _ len-1,而不是 len。
// 如果将其更改回 len,则将计算总体方差,而不是样本方差。
//
let denom = (self.len() - 1) as f64;
v / denom
}
}
fn std_dev(&self) -> f64 {
self.var().sqrt()
}
fn std_dev_pct(&self) -> f64 {
let hundred = 100_f64;
(self.std_dev() / self.mean()) * hundred
}
fn median_abs_dev(&self) -> f64 {
let med = self.median();
let abs_devs: Vec<f64> = self.iter().map(|&v| (med - v).abs()).collect();
// 这个常量是由比我更聪明的统计头脑得出的,但它与 R 和其他软件包对 MAD 的处理方式一致。
//
let number = 1.4826;
abs_devs.median() * number
}
fn median_abs_dev_pct(&self) -> f64 {
let hundred = 100_f64;
(self.median_abs_dev() / self.median()) * hundred
}
fn percentile(&self, pct: f64) -> f64 {
let mut tmp = self.to_vec();
local_sort(&mut tmp);
percentile_of_sorted(&tmp, pct)
}
fn quartiles(&self) -> (f64, f64, f64) {
let mut tmp = self.to_vec();
local_sort(&mut tmp);
let first = 25_f64;
let a = percentile_of_sorted(&tmp, first);
let second = 50_f64;
let b = percentile_of_sorted(&tmp, second);
let third = 75_f64;
let c = percentile_of_sorted(&tmp, third);
(a, b, c)
}
fn iqr(&self) -> f64 {
let (a, _, c) = self.quartiles();
c - a
}
}
// Helper 函数:使用线性插值提取代表排序后的样本集的 `pct` 百分位数的值。
// 如果样品未排序,则返回无意义的值。
fn percentile_of_sorted(sorted_samples: &[f64], pct: f64) -> f64 {
assert!(!sorted_samples.is_empty());
if sorted_samples.len() == 1 {
return sorted_samples[0];
}
let zero: f64 = 0.0;
assert!(zero <= pct);
let hundred = 100_f64;
assert!(pct <= hundred);
if pct == hundred {
return sorted_samples[sorted_samples.len() - 1];
}
let length = (sorted_samples.len() - 1) as f64;
let rank = (pct / hundred) * length;
let lrank = rank.floor();
let d = rank - lrank;
let n = lrank as usize;
let lo = sorted_samples[n];
let hi = sorted_samples[n + 1];
lo + (hi - lo) * d
}
/// Winsorize 一组样本,用这些百分数本身替换 `100-pct` 百分数以上和 `pct` 百分数以下的值。
/// 这是最小化离群值影响的方法,但要以使样本产生偏差为代价。
/// 它与修整的不同之处在于,它不会更改样本数,而只会更改异常值的值。
///
/// 请参见:<https://en.wikipedia.org/wiki/Winsorising>
///
///
pub fn winsorize(samples: &mut [f64], pct: f64) {
let mut tmp = samples.to_vec();
local_sort(&mut tmp);
let lo = percentile_of_sorted(&tmp, pct);
let hundred = 100_f64;
let hi = percentile_of_sorted(&tmp, hundred - pct);
for samp in samples {
if *samp > hi {
*samp = hi
} else if *samp < lo {
*samp = lo
}
}
}