Repository: gonum/gonum
Pull request #1570 (2021-02-20 15:43:57)

stat: add population (biased) variants of Variance, StdDev, MeanVariance and MeanStdDev
Please take a look.

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->


codecov-io: # [Codecov](https://codecov.io/gh/gonum/gonum/pull/1570?src=pr&el=h1) Report
> Merging [#1570](https://codecov.io/gh/gonum/gonum/pull/1570?src=pr&el=desc) (51819e2) into [master](https://codecov.io/gh/gonum/gonum/commit/3d9a0707517a8e08bfe9d2f3ea72ba7176bb84d6?el=desc) (3d9a070) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/gonum/gonum/pull/1570/graphs/tree.svg?width=650&height=150&src=pr&token=B7qHnqSPNL)](https://codecov.io/gh/gonum/gonum/pull/1570?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master    #1570   +/-   ##
=======================================
  Coverage   72.67%   72.68%           
=======================================
  Files         481      481           
  Lines       57919    57936   +17     
=======================================
+ Hits        42091    42108   +17     
  Misses      13428    13428           
  Partials     2400     2400           
```


| [Impacted Files](https://codecov.io/gh/gonum/gonum/pull/1570?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [stat/stat.go](https://codecov.io/gh/gonum/gonum/pull/1570/diff?src=pr&el=tree#diff-c3RhdC9zdGF0Lmdv) | `94.74% <100.00%> (+0.15%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gonum/gonum/pull/1570?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gonum/gonum/pull/1570?src=pr&el=footer). Last update [3d9a070...2eecd56](https://codecov.io/gh/gonum/gonum/pull/1570?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

kortschak: Good to see you back.
rilwen: I had to amend the last commit message. Can you reapprove?
romanwerpachowski: Good to be back! I had to amend the last commit message. Can you reapprove, please?
romanwerpachowski: Will do.

Patch:
diff --git a/AUTHORS b/AUTHORS
index 846c5f72..2429f44f 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -70,6 +70,7 @@ Mark Skilbeck <markskilbeck@gmail.com>
 Martin Diz <github@martindiz.com.ar>
 Matthieu Di Mercurio <matthieu.dimercurio@gmail.com>
 Max Halford <maxhalford25@gmail.com>
+Microsoft Corporation
 MinJae Kwon <k239507@gmail.com>
 Nathan Edwards <etaoinshrdluwho@gmail.com>
 Nick Potts <nick@the-potts.com>
diff --git a/stat/stat.go b/stat/stat.go
index e33188e5..bc38bb0b 100644
--- a/stat/stat.go
+++ b/stat/stat.go
@@ -1261,6 +1261,54 @@ func Variance(x, weights []float64) float64 {
 // len(x) must equal len(weights).
 // When weights sum to 1 or less, a biased variance estimator should be used.
 func MeanVariance(x, weights []float64) (mean, variance float64) {
+	var (
+		unnormalisedVariance float64
+		sumWeights           float64
+	)
+	mean, unnormalisedVariance, sumWeights = meanUnnormalisedVarianceSumWeights(x, weights)
+	return mean, unnormalisedVariance / (sumWeights - 1)
+}
+
+// PopMeanVariance computes the sample mean and biased variance (also known as
+// "population variance"), where the mean and variance are
+//  \sum_i w_i * x_i / (sum_i w_i)
+//  \sum_i w_i (x_i - mean)^2 / (sum_i w_i)
+// respectively.
+// If weights is nil then all of the weights are 1. If weights is not nil, then
+// len(x) must equal len(weights).
+func PopMeanVariance(x, weights []float64) (mean, variance float64) {
+	var (
+		unnormalisedVariance float64
+		sumWeights           float64
+	)
+	mean, unnormalisedVariance, sumWeights = meanUnnormalisedVarianceSumWeights(x, weights)
+	return mean, unnormalisedVariance / sumWeights
+}
+
+// PopMeanStdDev returns the sample mean and biased standard deviation
+// (also known as "population standard deviation").
+func PopMeanStdDev(x, weights []float64) (mean, std float64) {
+	mean, variance := PopMeanVariance(x, weights)
+	return mean, math.Sqrt(variance)
+}
+
+// PopStdDev returns the population standard deviation, i.e., a square root
+// of the biased variance estimate.
+func PopStdDev(x, weights []float64) float64 {
+	_, stDev := PopMeanStdDev(x, weights)
+	return stDev
+}
+
+// PopVariance computes the unbiased weighted sample variance:
+//  \sum_i w_i (x_i - mean)^2 / (sum_i w_i)
+// If weights is nil then all of the weights are 1. If weights is not nil, then
+// len(x) must equal len(weights).
+func PopVariance(x, weights []float64) float64 {
+	_, variance := PopMeanVariance(x, weights)
+	return variance
+}
+
+func meanUnnormalisedVarianceSumWeights(x, weights []float64) (mean, unnormalisedVariance, sumWeights float64) {
 	// This uses the corrected two-pass algorithm (1.7), from "Algorithms for computing
 	// the sample variance: Analysis and recommendations" by Chan, Tony F., Gene H. Golub,
 	// and Randall J. LeVeque.
@@ -1277,11 +1325,10 @@ func MeanVariance(x, weights []float64) (mean, variance float64) {
 			ss += d * d
 			compensation += d
 		}
-		variance = (ss - compensation*compensation/float64(len(x))) / float64(len(x)-1)
-		return mean, variance
+		unnormalisedVariance = (ss - compensation*compensation/float64(len(x)))
+		return mean, unnormalisedVariance, float64(len(x))
 	}
 
-	var sumWeights float64
 	for i, v := range x {
 		w := weights[i]
 		d := v - mean
@@ -1290,6 +1337,6 @@ func MeanVariance(x, weights []float64) (mean, variance float64) {
 		compensation += wd
 		sumWeights += w
 	}
-	variance = (ss - compensation*compensation/sumWeights) / (sumWeights - 1)
-	return mean, variance
+	unnormalisedVariance = (ss - compensation*compensation/sumWeights)
+	return mean, unnormalisedVariance, sumWeights
 }


Test patch:
diff --git a/stat/stat_test.go b/stat/stat_test.go
index c5cfa1a8..3b065753 100644
--- a/stat/stat_test.go
+++ b/stat/stat_test.go
@@ -1556,6 +1556,19 @@ func ExampleStdDev() {
 	// The weighted standard deviation of the samples is 10.5733
 }
 
+func ExamplePopStdDev() {
+	x := []float64{8, 2, -9, 15, 4}
+	stdev := PopStdDev(x, nil)
+	fmt.Printf("The standard deviation of the population is %.4f\n", stdev)
+
+	weights := []float64{2, 2, 6, 7, 1}
+	weightedStdev := PopStdDev(x, weights)
+	fmt.Printf("The weighted standard deviation of the population is %.4f\n", weightedStdev)
+	// Output:
+	// The standard deviation of the population is 7.8740
+	// The weighted standard deviation of the population is 10.2754
+}
+
 func ExampleStdErr() {
 	x := []float64{8, 2, -9, 15, 4}
 	weights := []float64{2, 2, 6, 7, 1}
@@ -1731,7 +1744,58 @@ func TestVariance(t *testing.T) {
 	if !panics(func() { Variance(make([]float64, 3), make([]float64, 2)) }) {
 		t.Errorf("Variance did not panic with x, weights length mismatch")
 	}
+}
 
+func TestPopVariance(t *testing.T) {
+	for i, test := range []struct {
+		x       []float64
+		weights []float64
+		ans     float64
+	}{
+		{
+			x:       []float64{8, -3, 7, 8, -4},
+			weights: nil,
+			ans:     30.16,
+		},
+		{
+			x:       []float64{8, -3, 7, 8, -4},
+			weights: []float64{1, 1, 1, 1, 1},
+			ans:     30.16,
+		},
+		{
+			x:       []float64{8, 3, 7, 8, 4},
+			weights: []float64{2, 1, 2, 1, 1},
+			ans:     3.6734693877551026,
+		},
+		{
+			x:       []float64{1, 4, 9},
+			weights: []float64{1, 1.5, 1},
+			ans:     9.387755102040817,
+		},
+		{
+			x:       []float64{1, 2, 3},
+			weights: []float64{1, 1.5, 1},
+			ans:     0.5714285714285714,
+		},
+		{
+			x:       []float64{2},
+			weights: nil,
+			ans:     0,
+		},
+		{
+			x:       []float64{2},
+			weights: []float64{2},
+			ans:     0,
+		},
+	} {
+		variance := PopVariance(test.x, test.weights)
+		if math.Abs(variance-test.ans) > 1e-14 {
+			t.Errorf("PopVariance mismatch case %d. Expected %v, Found %v", i, test.ans, variance)
+		}
+	}
+	if !panics(func() { PopVariance(make([]float64, 3), make([]float64, 2)) }) {
+		t.Errorf("PopVariance did not panic with x, weights length mismatch")
+	}
 }
 
 func ExampleVariance() {
@@ -1747,6 +1811,19 @@ func ExampleVariance() {
 	// The weighted variance of the samples is 111.7941
 }
 
+func ExamplePopVariance() {
+	x := []float64{8, 2, -9, 15, 4}
+	variance := PopVariance(x, nil)
+	fmt.Printf("The biased variance of the samples is %.4f\n", variance)
+
+	weights := []float64{2, 2, 6, 7, 1}
+	weightedVariance := PopVariance(x, weights)
+	fmt.Printf("The weighted biased variance of the samples is %.4f\n", weightedVariance)
+	// Output:
+	// The biased variance of the samples is 62.0000
+	// The weighted biased variance of the samples is 105.5833
+}
+
 func TestStdScore(t *testing.T) {
 	for i, test := range []struct {
 		x float64