Remove inferior estimate of error of weighted average
It was a biased estimator, while the alternative one isn't.
This commit is contained in:
parent
efaca98d0c
commit
77fa8b4ed2
@ -2,8 +2,8 @@
|
|||||||
//! sequence of numbers, and for their standard errors. The typical workflow
|
//! sequence of numbers, and for their standard errors. The typical workflow
|
||||||
//! looks like this:
|
//! looks like this:
|
||||||
//!
|
//!
|
||||||
//! 1. Initialize your estimator of choice (`Average`, `WeightedAverage` or
|
//! 1. Initialize your estimator of choice (`Average` or `WeightedAverage`) with
|
||||||
//! `WeightedAverage2`) with `new()`.
|
//! `new()`.
|
||||||
//! 2. Add some subset (called "samples") of the sequence of numbers (called
|
//! 2. Add some subset (called "samples") of the sequence of numbers (called
|
||||||
//! "population") for which you want to estimate the average, using `add()`
|
//! "population") for which you want to estimate the average, using `add()`
|
||||||
//! or `collect()`.
|
//! or `collect()`.
|
||||||
@ -32,8 +32,6 @@ extern crate conv;
|
|||||||
#[macro_use] mod macros;
|
#[macro_use] mod macros;
|
||||||
mod average;
|
mod average;
|
||||||
mod weighted_average;
|
mod weighted_average;
|
||||||
mod weighted_average2;
|
|
||||||
|
|
||||||
pub use average::Average;
|
pub use average::Average;
|
||||||
pub use weighted_average::WeightedAverage;
|
pub use weighted_average::WeightedAverage;
|
||||||
pub use weighted_average2::WeightedAverage as WeightedAverage2;
|
|
||||||
|
@ -1,7 +1,9 @@
|
|||||||
use core;
|
use core;
|
||||||
|
|
||||||
/// Estimate the weighted arithmetic mean and the weighted variance of a
|
use conv::ApproxFrom;
|
||||||
/// sequence of numbers ("population").
|
|
||||||
|
/// Estimate the weighted and unweighted arithmetic mean and the unweighted
|
||||||
|
/// variance of a sequence of numbers ("population").
|
||||||
///
|
///
|
||||||
/// This can be used to estimate the standard error of the weighted mean.
|
/// This can be used to estimate the standard error of the weighted mean.
|
||||||
///
|
///
|
||||||
@ -13,41 +15,58 @@ use core;
|
|||||||
///
|
///
|
||||||
/// let a: WeightedAverage = (1..6).zip(1..6)
|
/// let a: WeightedAverage = (1..6).zip(1..6)
|
||||||
/// .map(|(x, w)| (f64::from(x), f64::from(w))).collect();
|
/// .map(|(x, w)| (f64::from(x), f64::from(w))).collect();
|
||||||
/// println!("The weighted average is {} ± {}.", a.mean(), a.error());
|
/// println!("The weighted average is {} ± {}.", a.weighted_mean(), a.error());
|
||||||
/// ```
|
/// ```
|
||||||
#[derive(Debug, Clone)]
|
#[derive(Debug, Clone)]
|
||||||
pub struct WeightedAverage {
|
pub struct WeightedAverage {
|
||||||
/// Sum of the weights.
|
/// Sum of the weights.
|
||||||
weight_sum: f64,
|
weight_sum: f64,
|
||||||
/// Average value.
|
/// Sum of the squares of the weights.
|
||||||
avg: f64,
|
weight_sum_sq: f64,
|
||||||
/// Intermediate sum of squares for calculating the variance.
|
/// Weighted average value.
|
||||||
|
weighted_avg: f64,
|
||||||
|
|
||||||
|
/// Number of samples.
|
||||||
|
n: u64,
|
||||||
|
/// Unweighted average value.
|
||||||
|
unweighted_avg: f64,
|
||||||
|
/// Intermediate sum of squares for calculating the *unweighted* variance.
|
||||||
v: f64,
|
v: f64,
|
||||||
}
|
}
|
||||||
|
|
||||||
impl WeightedAverage {
|
impl WeightedAverage {
|
||||||
/// Create a new weighted average estimator.
|
/// Create a new weighted and unweighted average estimator.
|
||||||
pub fn new() -> WeightedAverage {
|
pub fn new() -> WeightedAverage {
|
||||||
WeightedAverage { weight_sum: 0., avg: 0., v: 0. }
|
WeightedAverage {
|
||||||
|
weight_sum: 0., weight_sum_sq: 0., weighted_avg: 0.,
|
||||||
|
n: 0, unweighted_avg: 0., v: 0.,
|
||||||
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Add a weighted element sampled from the population.
|
/// Add a weighted element sampled from the population.
|
||||||
pub fn add(&mut self, sample: f64, weight: f64) {
|
pub fn add(&mut self, sample: f64, weight: f64) {
|
||||||
// This algorithm was suggested by West in 1979.
|
// The algorithm for the unweighted average was suggested by Welford in 1962.
|
||||||
|
// The algorithm for the weighted average was suggested by West in 1979.
|
||||||
//
|
//
|
||||||
// See
|
// See
|
||||||
// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
|
// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
|
||||||
// and
|
// and
|
||||||
// http://people.ds.cam.ac.uk/fanf2/hermes/doc/antiforgery/stats.pdf.
|
// http://people.ds.cam.ac.uk/fanf2/hermes/doc/antiforgery/stats.pdf.
|
||||||
self.weight_sum += weight;
|
self.weight_sum += weight;
|
||||||
let prev_avg = self.avg;
|
self.weight_sum_sq += weight*weight;
|
||||||
self.avg = prev_avg + (weight / self.weight_sum) * (sample - prev_avg);
|
|
||||||
self.v += weight * (sample - prev_avg) * (sample - self.avg);
|
let prev_avg = self.weighted_avg;
|
||||||
|
self.weighted_avg = prev_avg + (weight / self.weight_sum) * (sample - prev_avg);
|
||||||
|
|
||||||
|
self.n += 1;
|
||||||
|
let delta = sample - self.unweighted_avg;
|
||||||
|
self.unweighted_avg += delta / f64::approx_from(self.n).unwrap();
|
||||||
|
self.v += delta * (sample - self.unweighted_avg);
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Determine whether the samples are empty.
|
/// Determine whether the sample is empty.
|
||||||
pub fn is_empty(&self) -> bool {
|
pub fn is_empty(&self) -> bool {
|
||||||
self.weight_sum == 0. && self.v == 0. && self.avg == 0.
|
self.n == 0
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Return the sum of the weights.
|
/// Return the sum of the weights.
|
||||||
@ -55,57 +74,70 @@ impl WeightedAverage {
|
|||||||
self.weight_sum
|
self.weight_sum
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Estimate the weighted mean of the population.
|
/// Return the sum of the squared weights.
|
||||||
pub fn mean(&self) -> f64 {
|
pub fn sum_weights_sq(&self) -> f64 {
|
||||||
self.avg
|
self.weight_sum_sq
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Calculate the weighted population variance of the sample.
|
/// Estimate the weighted mean of the sequence.
|
||||||
///
|
pub fn weighted_mean(&self) -> f64 {
|
||||||
/// This is a biased estimator of the weighted variance of the population.
|
self.weighted_avg
|
||||||
pub fn population_variance(&self) -> f64 {
|
}
|
||||||
|
|
||||||
|
/// Estimate the unweighted mean of the sequence.
|
||||||
|
pub fn unweighted_mean(&self) -> f64 {
|
||||||
|
self.unweighted_avg
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Return sample size.
|
||||||
|
pub fn len(&self) -> u64 {
|
||||||
|
self.n
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Calculate the effective sample size.
|
||||||
|
pub fn effective_len(&self) -> f64 {
|
||||||
if self.is_empty() {
|
if self.is_empty() {
|
||||||
0.
|
return 0.
|
||||||
} else {
|
|
||||||
self.v / self.weight_sum
|
|
||||||
}
|
}
|
||||||
|
self.weight_sum * self.weight_sum / self.weight_sum_sq
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Calculate the weighted sample variance.
|
/// Calculate the *unweighted* population variance of the sample.
|
||||||
///
|
///
|
||||||
/// This is an unbiased estimator of the weighted variance of the population.
|
/// This is a biased estimator of the variance of the population.
|
||||||
|
pub fn population_variance(&self) -> f64 {
|
||||||
|
if self.n < 2 {
|
||||||
|
return 0.;
|
||||||
|
}
|
||||||
|
self.v / f64::approx_from(self.n).unwrap()
|
||||||
|
}
|
||||||
|
|
||||||
|
/// Calculate the *unweighted* sample variance.
|
||||||
///
|
///
|
||||||
/// Note that this will return 0 if the sum of the weights is <= 1.
|
/// This is an unbiased estimator of the variance of the population.
|
||||||
pub fn sample_variance(&self) -> f64 {
|
pub fn sample_variance(&self) -> f64 {
|
||||||
if self.weight_sum <= 1. {
|
if self.n < 2 {
|
||||||
0.
|
return 0.;
|
||||||
} else {
|
|
||||||
self.v / (self.weight_sum - 1.0)
|
|
||||||
}
|
}
|
||||||
|
self.v / f64::approx_from(self.n - 1).unwrap()
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Estimate the standard error of the weighted mean of the population.
|
/// Estimate the standard error of the *weighted* mean of the sequence.
|
||||||
///
|
///
|
||||||
/// Note that this will return 0 if the sum of the weights is 0.
|
/// Returns 0 if the sum of weights is 0.
|
||||||
/// For this estimator, the sum of weights should be larger than 1.
|
|
||||||
///
|
///
|
||||||
/// This biased estimator uses the weighted variance and the sum of weights.
|
/// This unbiased estimator assumes that the samples were independently
|
||||||
/// It considers the weights as (noninteger) counts of how often the sample
|
/// drawn from the same population with constant variance.
|
||||||
/// has been observed, applying the standard formulas to calculate mean,
|
|
||||||
/// variance and sample size across all "repeats".
|
|
||||||
pub fn error(&self) -> f64 {
|
pub fn error(&self) -> f64 {
|
||||||
// This uses the same estimate as SPSS.
|
// This uses the same estimate as WinCross, which should provide better
|
||||||
|
// results than the ones used by SPSS or Mentor.
|
||||||
//
|
//
|
||||||
// See http://www.analyticalgroup.com/download/WEIGHTED_MEAN.pdf.
|
// See http://www.analyticalgroup.com/download/WEIGHTED_VARIANCE.pdf.
|
||||||
if self.weight_sum == 0. {
|
if self.weight_sum == 0. {
|
||||||
return 0.;
|
return 0.;
|
||||||
}
|
}
|
||||||
let variance = if self.weight_sum <= 1. {
|
let inv_effective_len = self.weight_sum_sq / (self.weight_sum * self.weight_sum);
|
||||||
self.population_variance()
|
(self.sample_variance() * inv_effective_len).sqrt()
|
||||||
} else {
|
|
||||||
self.sample_variance()
|
|
||||||
};
|
|
||||||
(variance / self.weight_sum).sqrt()
|
|
||||||
}
|
}
|
||||||
|
|
||||||
/// Merge another sample into this one.
|
/// Merge another sample into this one.
|
||||||
@ -124,20 +156,32 @@ impl WeightedAverage {
|
|||||||
/// let mut avg_left: WeightedAverage = left.iter().map(|&x| x).collect();
|
/// let mut avg_left: WeightedAverage = left.iter().map(|&x| x).collect();
|
||||||
/// let avg_right: WeightedAverage = right.iter().map(|&x| x).collect();
|
/// let avg_right: WeightedAverage = right.iter().map(|&x| x).collect();
|
||||||
/// avg_left.merge(&avg_right);
|
/// avg_left.merge(&avg_right);
|
||||||
/// assert!((avg_total.mean() - avg_left.mean()).abs() < 1e-15);
|
/// assert!((avg_total.weighted_mean() - avg_left.weighted_mean()).abs() < 1e-15);
|
||||||
/// assert!((avg_total.error() - avg_left.error()).abs() < 1e-15);
|
/// assert!((avg_total.error() - avg_left.error()).abs() < 1e-15);
|
||||||
/// ```
|
/// ```
|
||||||
pub fn merge(&mut self, other: &WeightedAverage) {
|
pub fn merge(&mut self, other: &WeightedAverage) {
|
||||||
// This is similar to the algorithm proposed by Chan et al. in 1979.
|
// This is similar to the algorithm proposed by Chan et al. in 1979.
|
||||||
//
|
//
|
||||||
// See https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance.
|
// See https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance.
|
||||||
let delta = other.avg - self.avg;
|
{
|
||||||
let total_weight_sum = self.weight_sum + other.weight_sum;
|
let total_weight_sum = self.weight_sum + other.weight_sum;
|
||||||
self.avg = (self.weight_sum * self.avg + other.weight_sum * other.avg)
|
self.weighted_avg = (self.weight_sum * self.weighted_avg
|
||||||
|
+ other.weight_sum * other.weighted_avg)
|
||||||
/ (self.weight_sum + other.weight_sum);
|
/ (self.weight_sum + other.weight_sum);
|
||||||
self.v += other.v + delta*delta * self.weight_sum * other.weight_sum
|
|
||||||
/ total_weight_sum;
|
|
||||||
self.weight_sum = total_weight_sum;
|
self.weight_sum = total_weight_sum;
|
||||||
|
self.weight_sum_sq += other.weight_sum_sq;
|
||||||
|
}
|
||||||
|
{
|
||||||
|
let delta = other.unweighted_avg - self.unweighted_avg;
|
||||||
|
let len_self = f64::approx_from(self.n).unwrap();
|
||||||
|
let len_other = f64::approx_from(other.n).unwrap();
|
||||||
|
let len_total = len_self + len_other;
|
||||||
|
self.n += other.n;
|
||||||
|
self.unweighted_avg = (len_self * self.unweighted_avg
|
||||||
|
+ len_other * other.unweighted_avg)
|
||||||
|
/ len_total;
|
||||||
|
self.v += other.v + delta*delta * len_self * len_other / len_total;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -172,8 +216,11 @@ mod tests {
|
|||||||
let mut avg_left: WeightedAverage = left.iter().map(|x| (*x, 1.)).collect();
|
let mut avg_left: WeightedAverage = left.iter().map(|x| (*x, 1.)).collect();
|
||||||
let avg_right: WeightedAverage = right.iter().map(|x| (*x, 1.)).collect();
|
let avg_right: WeightedAverage = right.iter().map(|x| (*x, 1.)).collect();
|
||||||
avg_left.merge(&avg_right);
|
avg_left.merge(&avg_right);
|
||||||
|
assert_eq!(avg_total.n, avg_left.n);
|
||||||
assert_eq!(avg_total.weight_sum, avg_left.weight_sum);
|
assert_eq!(avg_total.weight_sum, avg_left.weight_sum);
|
||||||
assert_eq!(avg_total.avg, avg_left.avg);
|
assert_eq!(avg_total.weight_sum_sq, avg_left.weight_sum_sq);
|
||||||
|
assert_eq!(avg_total.weighted_avg, avg_left.weighted_avg);
|
||||||
|
assert_eq!(avg_total.unweighted_avg, avg_left.unweighted_avg);
|
||||||
assert_eq!(avg_total.v, avg_left.v);
|
assert_eq!(avg_total.v, avg_left.v);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
@ -189,8 +236,11 @@ mod tests {
|
|||||||
let mut avg_left: WeightedAverage = left.iter().map(|&(x, w)| (x, w)).collect();
|
let mut avg_left: WeightedAverage = left.iter().map(|&(x, w)| (x, w)).collect();
|
||||||
let avg_right: WeightedAverage = right.iter().map(|&(x, w)| (x, w)).collect();
|
let avg_right: WeightedAverage = right.iter().map(|&(x, w)| (x, w)).collect();
|
||||||
avg_left.merge(&avg_right);
|
avg_left.merge(&avg_right);
|
||||||
|
assert_eq!(avg_total.n, avg_left.n);
|
||||||
assert_almost_eq!(avg_total.weight_sum, avg_left.weight_sum, 1e-15);
|
assert_almost_eq!(avg_total.weight_sum, avg_left.weight_sum, 1e-15);
|
||||||
assert_almost_eq!(avg_total.avg, avg_left.avg, 1e-15);
|
assert_eq!(avg_total.weight_sum_sq, avg_left.weight_sum_sq);
|
||||||
|
assert_almost_eq!(avg_total.weighted_avg, avg_left.weighted_avg, 1e-15);
|
||||||
|
assert_almost_eq!(avg_total.unweighted_avg, avg_left.unweighted_avg, 1e-15);
|
||||||
assert_almost_eq!(avg_total.v, avg_left.v, 1e-14);
|
assert_almost_eq!(avg_total.v, avg_left.v, 1e-14);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -1,246 +0,0 @@
|
|||||||
use core;
|
|
||||||
|
|
||||||
use conv::ApproxFrom;
|
|
||||||
|
|
||||||
/// Estimate the weighted and unweighted arithmetic mean and the unweighted
|
|
||||||
/// variance of a sequence of numbers ("population").
|
|
||||||
///
|
|
||||||
/// This can be used to estimate the standard error of the weighted mean.
|
|
||||||
///
|
|
||||||
///
|
|
||||||
/// ## Example
|
|
||||||
///
|
|
||||||
/// ```
|
|
||||||
/// use average::WeightedAverage2 as WeightedAverage;
|
|
||||||
///
|
|
||||||
/// let a: WeightedAverage = (1..6).zip(1..6)
|
|
||||||
/// .map(|(x, w)| (f64::from(x), f64::from(w))).collect();
|
|
||||||
/// println!("The weighted average is {} ± {}.", a.weighted_mean(), a.error());
|
|
||||||
/// ```
|
|
||||||
#[derive(Debug, Clone)]
|
|
||||||
pub struct WeightedAverage {
|
|
||||||
/// Sum of the weights.
|
|
||||||
weight_sum: f64,
|
|
||||||
/// Sum of the squares of the weights.
|
|
||||||
weight_sum_sq: f64,
|
|
||||||
/// Weighted average value.
|
|
||||||
weighted_avg: f64,
|
|
||||||
|
|
||||||
/// Number of samples.
|
|
||||||
n: u64,
|
|
||||||
/// Unweighted average value.
|
|
||||||
unweighted_avg: f64,
|
|
||||||
/// Intermediate sum of squares for calculating the *unweighted* variance.
|
|
||||||
v: f64,
|
|
||||||
}
|
|
||||||
|
|
||||||
impl WeightedAverage {
|
|
||||||
/// Create a new weighted and unweighted average estimator.
|
|
||||||
pub fn new() -> WeightedAverage {
|
|
||||||
WeightedAverage {
|
|
||||||
weight_sum: 0., weight_sum_sq: 0., weighted_avg: 0.,
|
|
||||||
n: 0, unweighted_avg: 0., v: 0.,
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Add a weighted element sampled from the population.
|
|
||||||
pub fn add(&mut self, sample: f64, weight: f64) {
|
|
||||||
// The algorithm for the unweighted average was suggested by Welford in 1962.
|
|
||||||
// The algorithm for the weighted average was suggested by West in 1979.
|
|
||||||
//
|
|
||||||
// See
|
|
||||||
// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
|
|
||||||
// and
|
|
||||||
// http://people.ds.cam.ac.uk/fanf2/hermes/doc/antiforgery/stats.pdf.
|
|
||||||
self.weight_sum += weight;
|
|
||||||
self.weight_sum_sq += weight*weight;
|
|
||||||
|
|
||||||
let prev_avg = self.weighted_avg;
|
|
||||||
self.weighted_avg = prev_avg + (weight / self.weight_sum) * (sample - prev_avg);
|
|
||||||
|
|
||||||
self.n += 1;
|
|
||||||
let delta = sample - self.unweighted_avg;
|
|
||||||
self.unweighted_avg += delta / f64::approx_from(self.n).unwrap();
|
|
||||||
self.v += delta * (sample - self.unweighted_avg);
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Determine whether the sample is empty.
|
|
||||||
pub fn is_empty(&self) -> bool {
|
|
||||||
self.n == 0
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Return the sum of the weights.
|
|
||||||
pub fn sum_weights(&self) -> f64 {
|
|
||||||
self.weight_sum
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Return the sum of the squared weights.
|
|
||||||
pub fn sum_weights_sq(&self) -> f64 {
|
|
||||||
self.weight_sum_sq
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Estimate the weighted mean of the sequence.
|
|
||||||
pub fn weighted_mean(&self) -> f64 {
|
|
||||||
self.weighted_avg
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Estimate the unweighted mean of the sequence.
|
|
||||||
pub fn unweighted_mean(&self) -> f64 {
|
|
||||||
self.unweighted_avg
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Return sample size.
|
|
||||||
pub fn len(&self) -> u64 {
|
|
||||||
self.n
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Calculate the effective sample size.
|
|
||||||
pub fn effective_len(&self) -> f64 {
|
|
||||||
if self.is_empty() {
|
|
||||||
return 0.
|
|
||||||
}
|
|
||||||
self.weight_sum * self.weight_sum / self.weight_sum_sq
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Calculate the *unweighted* population variance of the sample.
|
|
||||||
///
|
|
||||||
/// This is a biased estimator of the variance of the population.
|
|
||||||
pub fn population_variance(&self) -> f64 {
|
|
||||||
if self.n < 2 {
|
|
||||||
return 0.;
|
|
||||||
}
|
|
||||||
self.v / f64::approx_from(self.n).unwrap()
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Calculate the *unweighted* sample variance.
|
|
||||||
///
|
|
||||||
/// This is an unbiased estimator of the variance of the population.
|
|
||||||
pub fn sample_variance(&self) -> f64 {
|
|
||||||
if self.n < 2 {
|
|
||||||
return 0.;
|
|
||||||
}
|
|
||||||
self.v / f64::approx_from(self.n - 1).unwrap()
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Estimate the standard error of the *weighted* mean of the sequence.
|
|
||||||
///
|
|
||||||
/// Returns 0 if the sum of weights is 0.
|
|
||||||
///
|
|
||||||
/// This unbiased estimator assumes that the samples were independently
|
|
||||||
/// drawn from the same population with constant variance.
|
|
||||||
pub fn error(&self) -> f64 {
|
|
||||||
// This uses the same estimate as WinCross.
|
|
||||||
//
|
|
||||||
// See http://www.analyticalgroup.com/download/WEIGHTED_MEAN.pdf.
|
|
||||||
if self.weight_sum_sq == 0. || self.weight_sum == 0. {
|
|
||||||
return 0.;
|
|
||||||
}
|
|
||||||
let effective_base = self.weight_sum * self.weight_sum / self.weight_sum_sq;
|
|
||||||
(self.sample_variance() / effective_base).sqrt()
|
|
||||||
}
|
|
||||||
|
|
||||||
/// Merge another sample into this one.
|
|
||||||
///
|
|
||||||
///
|
|
||||||
/// ## Example
|
|
||||||
///
|
|
||||||
/// ```
|
|
||||||
/// use average::WeightedAverage2 as WeightedAverage;
|
|
||||||
///
|
|
||||||
/// let weighted_sequence: &[(f64, f64)] = &[
|
|
||||||
/// (1., 0.1), (2., 0.2), (3., 0.3), (4., 0.4), (5., 0.5),
|
|
||||||
/// (6., 0.6), (7., 0.7), (8., 0.8), (9., 0.9)];
|
|
||||||
/// let (left, right) = weighted_sequence.split_at(3);
|
|
||||||
/// let avg_total: WeightedAverage = weighted_sequence.iter().map(|&x| x).collect();
|
|
||||||
/// let mut avg_left: WeightedAverage = left.iter().map(|&x| x).collect();
|
|
||||||
/// let avg_right: WeightedAverage = right.iter().map(|&x| x).collect();
|
|
||||||
/// avg_left.merge(&avg_right);
|
|
||||||
/// assert!((avg_total.weighted_mean() - avg_left.weighted_mean()).abs() < 1e-15);
|
|
||||||
/// assert!((avg_total.error() - avg_left.error()).abs() < 1e-15);
|
|
||||||
/// ```
|
|
||||||
pub fn merge(&mut self, other: &WeightedAverage) {
|
|
||||||
// This is similar to the algorithm proposed by Chan et al. in 1979.
|
|
||||||
//
|
|
||||||
// See https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance.
|
|
||||||
{
|
|
||||||
let total_weight_sum = self.weight_sum + other.weight_sum;
|
|
||||||
self.weighted_avg = (self.weight_sum * self.weighted_avg
|
|
||||||
+ other.weight_sum * other.weighted_avg)
|
|
||||||
/ (self.weight_sum + other.weight_sum);
|
|
||||||
self.weight_sum = total_weight_sum;
|
|
||||||
self.weight_sum_sq += other.weight_sum_sq;
|
|
||||||
}
|
|
||||||
{
|
|
||||||
let delta = other.unweighted_avg - self.unweighted_avg;
|
|
||||||
let len_self = f64::approx_from(self.n).unwrap();
|
|
||||||
let len_other = f64::approx_from(other.n).unwrap();
|
|
||||||
let len_total = len_self + len_other;
|
|
||||||
self.n += other.n;
|
|
||||||
self.unweighted_avg = (len_self * self.unweighted_avg
|
|
||||||
+ len_other * other.unweighted_avg)
|
|
||||||
/ len_total;
|
|
||||||
self.v += other.v + delta*delta * len_self * len_other / len_total;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
impl core::default::Default for WeightedAverage {
|
|
||||||
fn default() -> WeightedAverage {
|
|
||||||
WeightedAverage::new()
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
impl core::iter::FromIterator<(f64, f64)> for WeightedAverage {
|
|
||||||
fn from_iter<T>(iter: T) -> WeightedAverage
|
|
||||||
where T: IntoIterator<Item=(f64, f64)>
|
|
||||||
{
|
|
||||||
let mut a = WeightedAverage::new();
|
|
||||||
for (i, w) in iter {
|
|
||||||
a.add(i, w);
|
|
||||||
}
|
|
||||||
a
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
#[cfg(test)]
|
|
||||||
mod tests {
|
|
||||||
use super::*;
|
|
||||||
|
|
||||||
#[test]
|
|
||||||
fn merge_unweighted() {
|
|
||||||
let sequence: &[f64] = &[1., 2., 3., 4., 5., 6., 7., 8., 9.];
|
|
||||||
for mid in 0..sequence.len() {
|
|
||||||
let (left, right) = sequence.split_at(mid);
|
|
||||||
let avg_total: WeightedAverage = sequence.iter().map(|x| (*x, 1.)).collect();
|
|
||||||
let mut avg_left: WeightedAverage = left.iter().map(|x| (*x, 1.)).collect();
|
|
||||||
let avg_right: WeightedAverage = right.iter().map(|x| (*x, 1.)).collect();
|
|
||||||
avg_left.merge(&avg_right);
|
|
||||||
assert_eq!(avg_total.n, avg_left.n);
|
|
||||||
assert_eq!(avg_total.weight_sum, avg_left.weight_sum);
|
|
||||||
assert_eq!(avg_total.weight_sum_sq, avg_left.weight_sum_sq);
|
|
||||||
assert_eq!(avg_total.weighted_avg, avg_left.weighted_avg);
|
|
||||||
assert_eq!(avg_total.unweighted_avg, avg_left.unweighted_avg);
|
|
||||||
assert_eq!(avg_total.v, avg_left.v);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
#[test]
|
|
||||||
fn merge_weighted() {
|
|
||||||
let sequence: &[(f64, f64)] = &[
|
|
||||||
(1., 0.1), (2., 0.2), (3., 0.3), (4., 0.4), (5., 0.5),
|
|
||||||
(6., 0.6), (7., 0.7), (8., 0.8), (9., 0.)];
|
|
||||||
for mid in 0..sequence.len() {
|
|
||||||
let (left, right) = sequence.split_at(mid);
|
|
||||||
let avg_total: WeightedAverage = sequence.iter().map(|&(x, w)| (x, w)).collect();
|
|
||||||
let mut avg_left: WeightedAverage = left.iter().map(|&(x, w)| (x, w)).collect();
|
|
||||||
let avg_right: WeightedAverage = right.iter().map(|&(x, w)| (x, w)).collect();
|
|
||||||
avg_left.merge(&avg_right);
|
|
||||||
assert_eq!(avg_total.n, avg_left.n);
|
|
||||||
assert_almost_eq!(avg_total.weight_sum, avg_left.weight_sum, 1e-15);
|
|
||||||
assert_eq!(avg_total.weight_sum_sq, avg_left.weight_sum_sq);
|
|
||||||
assert_almost_eq!(avg_total.weighted_avg, avg_left.weighted_avg, 1e-15);
|
|
||||||
assert_almost_eq!(avg_total.unweighted_avg, avg_left.unweighted_avg, 1e-15);
|
|
||||||
assert_almost_eq!(avg_total.v, avg_left.v, 1e-14);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
@ -9,15 +9,23 @@ use average::WeightedAverage;
|
|||||||
#[test]
|
#[test]
|
||||||
fn trivial() {
|
fn trivial() {
|
||||||
let mut a = WeightedAverage::new();
|
let mut a = WeightedAverage::new();
|
||||||
|
assert_eq!(a.len(), 0);
|
||||||
assert_eq!(a.sum_weights(), 0.);
|
assert_eq!(a.sum_weights(), 0.);
|
||||||
|
assert_eq!(a.sum_weights_sq(), 0.);
|
||||||
a.add(1.0, 1.0);
|
a.add(1.0, 1.0);
|
||||||
assert_eq!(a.mean(), 1.0);
|
assert_eq!(a.len(), 1);
|
||||||
|
assert_eq!(a.weighted_mean(), 1.0);
|
||||||
|
assert_eq!(a.unweighted_mean(), 1.0);
|
||||||
assert_eq!(a.sum_weights(), 1.0);
|
assert_eq!(a.sum_weights(), 1.0);
|
||||||
|
assert_eq!(a.sum_weights_sq(), 1.0);
|
||||||
assert_eq!(a.population_variance(), 0.0);
|
assert_eq!(a.population_variance(), 0.0);
|
||||||
assert_eq!(a.error(), 0.0);
|
assert_eq!(a.error(), 0.0);
|
||||||
a.add(1.0, 1.0);
|
a.add(1.0, 1.0);
|
||||||
assert_eq!(a.mean(), 1.0);
|
assert_eq!(a.len(), 2);
|
||||||
|
assert_eq!(a.weighted_mean(), 1.0);
|
||||||
|
assert_eq!(a.unweighted_mean(), 1.0);
|
||||||
assert_eq!(a.sum_weights(), 2.0);
|
assert_eq!(a.sum_weights(), 2.0);
|
||||||
|
assert_eq!(a.sum_weights_sq(), 2.0);
|
||||||
assert_eq!(a.population_variance(), 0.0);
|
assert_eq!(a.population_variance(), 0.0);
|
||||||
assert_eq!(a.error(), 0.0);
|
assert_eq!(a.error(), 0.0);
|
||||||
}
|
}
|
||||||
@ -25,7 +33,9 @@ fn trivial() {
|
|||||||
#[test]
|
#[test]
|
||||||
fn simple() {
|
fn simple() {
|
||||||
let a: WeightedAverage = (1..6).map(|x| (f64::from(x), 1.0)).collect();
|
let a: WeightedAverage = (1..6).map(|x| (f64::from(x), 1.0)).collect();
|
||||||
assert_eq!(a.mean(), 3.0);
|
assert_eq!(a.len(), 5);
|
||||||
|
assert_eq!(a.weighted_mean(), 3.0);
|
||||||
|
assert_eq!(a.unweighted_mean(), 3.0);
|
||||||
assert_eq!(a.sum_weights(), 5.0);
|
assert_eq!(a.sum_weights(), 5.0);
|
||||||
assert_eq!(a.sample_variance(), 2.5);
|
assert_eq!(a.sample_variance(), 2.5);
|
||||||
assert_almost_eq!(a.error(), f64::sqrt(0.5), 1e-16);
|
assert_almost_eq!(a.error(), f64::sqrt(0.5), 1e-16);
|
||||||
@ -38,10 +48,12 @@ fn reference() {
|
|||||||
let weights = &[1.23, 2.12, 1.23, 0.32, 1.53, 0.59, 0.94, 0.94, 0.84, 0.73];
|
let weights = &[1.23, 2.12, 1.23, 0.32, 1.53, 0.59, 0.94, 0.94, 0.84, 0.73];
|
||||||
let a: WeightedAverage = values.iter().zip(weights.iter())
|
let a: WeightedAverage = values.iter().zip(weights.iter())
|
||||||
.map(|(x, w)| (*x, *w)).collect();
|
.map(|(x, w)| (*x, *w)).collect();
|
||||||
assert_almost_eq!(a.mean(), 3.53486, 1e-5);
|
assert_almost_eq!(a.weighted_mean(), 3.53486, 1e-5);
|
||||||
assert_almost_eq!(a.sample_variance(), 1.8210, 1e-4);
|
assert_almost_eq!(a.sample_variance(), 1.7889, 1e-4);
|
||||||
assert_eq!(a.sum_weights(), 10.47);
|
assert_eq!(a.sum_weights(), 10.47);
|
||||||
assert_almost_eq!(a.error(), f64::sqrt(0.1739), 1e-4);
|
assert_eq!(a.len(), 10);
|
||||||
|
assert_almost_eq!(a.effective_len(), 8.2315, 1e-4);
|
||||||
|
assert_almost_eq!(a.error(), f64::sqrt(0.2173), 1e-4);
|
||||||
}
|
}
|
||||||
|
|
||||||
#[test]
|
#[test]
|
||||||
|
@ -1,66 +0,0 @@
|
|||||||
#[macro_use] extern crate average;
|
|
||||||
|
|
||||||
extern crate core;
|
|
||||||
|
|
||||||
use core::iter::Iterator;
|
|
||||||
|
|
||||||
use average::WeightedAverage2 as WeightedAverage;
|
|
||||||
|
|
||||||
#[test]
|
|
||||||
fn trivial() {
|
|
||||||
let mut a = WeightedAverage::new();
|
|
||||||
assert_eq!(a.len(), 0);
|
|
||||||
assert_eq!(a.sum_weights(), 0.);
|
|
||||||
assert_eq!(a.sum_weights_sq(), 0.);
|
|
||||||
a.add(1.0, 1.0);
|
|
||||||
assert_eq!(a.len(), 1);
|
|
||||||
assert_eq!(a.weighted_mean(), 1.0);
|
|
||||||
assert_eq!(a.unweighted_mean(), 1.0);
|
|
||||||
assert_eq!(a.sum_weights(), 1.0);
|
|
||||||
assert_eq!(a.sum_weights_sq(), 1.0);
|
|
||||||
assert_eq!(a.population_variance(), 0.0);
|
|
||||||
assert_eq!(a.error(), 0.0);
|
|
||||||
a.add(1.0, 1.0);
|
|
||||||
assert_eq!(a.len(), 2);
|
|
||||||
assert_eq!(a.weighted_mean(), 1.0);
|
|
||||||
assert_eq!(a.unweighted_mean(), 1.0);
|
|
||||||
assert_eq!(a.sum_weights(), 2.0);
|
|
||||||
assert_eq!(a.sum_weights_sq(), 2.0);
|
|
||||||
assert_eq!(a.population_variance(), 0.0);
|
|
||||||
assert_eq!(a.error(), 0.0);
|
|
||||||
}
|
|
||||||
|
|
||||||
#[test]
|
|
||||||
fn simple() {
|
|
||||||
let a: WeightedAverage = (1..6).map(|x| (f64::from(x), 1.0)).collect();
|
|
||||||
assert_eq!(a.len(), 5);
|
|
||||||
assert_eq!(a.weighted_mean(), 3.0);
|
|
||||||
assert_eq!(a.unweighted_mean(), 3.0);
|
|
||||||
assert_eq!(a.sum_weights(), 5.0);
|
|
||||||
assert_eq!(a.sample_variance(), 2.5);
|
|
||||||
assert_almost_eq!(a.error(), f64::sqrt(0.5), 1e-16);
|
|
||||||
}
|
|
||||||
|
|
||||||
#[test]
|
|
||||||
fn reference() {
|
|
||||||
// Example from http://www.analyticalgroup.com/download/WEIGHTED_MEAN.pdf.
|
|
||||||
let values = &[5., 5., 4., 4., 3., 4., 3., 2., 2., 1.];
|
|
||||||
let weights = &[1.23, 2.12, 1.23, 0.32, 1.53, 0.59, 0.94, 0.94, 0.84, 0.73];
|
|
||||||
let a: WeightedAverage = values.iter().zip(weights.iter())
|
|
||||||
.map(|(x, w)| (*x, *w)).collect();
|
|
||||||
assert_almost_eq!(a.weighted_mean(), 3.53486, 1e-5);
|
|
||||||
assert_almost_eq!(a.sample_variance(), 1.7889, 1e-4);
|
|
||||||
assert_eq!(a.sum_weights(), 10.47);
|
|
||||||
assert_eq!(a.len(), 10);
|
|
||||||
assert_almost_eq!(a.effective_len(), 8.2315, 1e-4);
|
|
||||||
assert_almost_eq!(a.error(), f64::sqrt(0.2173), 1e-4);
|
|
||||||
}
|
|
||||||
|
|
||||||
#[test]
|
|
||||||
fn error_corner_case() {
|
|
||||||
let values = &[1., 2.];
|
|
||||||
let weights = &[0.5, 0.5];
|
|
||||||
let a: WeightedAverage = values.iter().zip(weights.iter())
|
|
||||||
.map(|(x, w)| (*x, *w)).collect();
|
|
||||||
assert_eq!(a.error(), 0.5);
|
|
||||||
}
|
|
Loading…
Reference in New Issue
Block a user