diff --git a/covariancematrix.go b/covariancematrix.go new file mode 100644 index 00000000..06515ed6 --- /dev/null +++ b/covariancematrix.go @@ -0,0 +1,73 @@ +// Copyright ©2014 The gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package stat + +import ( + "github.com/gonum/floats" + "github.com/gonum/matrix/mat64" +) + +// CovarianceMatrix calculates a covariance matrix (also known as a +// variance-covariance matrix) from a matrix of data, using a two-pass +// algorithm. The matrix returned will be symmetric, square, and +// positive-semidefinite. +// +// The weights wts should have the length equal to the number of rows in +// input data matrix x. cov should either be a square matrix with the same +// number of columns as the input data matrix x, or nil in which case a new +// Dense matrix will be constructed. +func CovarianceMatrix(cov *mat64.Dense, x mat64.Matrix, wts []float64) *mat64.Dense { + // This is the matrix version of the two-pass algorithm. It doesn't use the + // additional floating point error correction that the Covariance function uses + // to reduce the impact of rounding during centering. + + r, c := x.Dims() + + // TODO(jonlawlor): indicate that the resulting matrix is symmetric, which + // should improve performance. + if cov == nil { + cov = mat64.NewDense(c, c, nil) + } else if covr, covc := cov.Dims(); covr != covc || covc != c { + panic(mat64.ErrShape) + } + + var xt mat64.Dense + xt.TCopy(x) + // Subtract the mean of each of the columns. + for i := 0; i < c; i++ { + v := xt.RowView(i) + mean := Mean(v, wts) + floats.AddConst(-mean, v) + } + + var xc mat64.Dense + xc.TCopy(&xt) + + var n, scale float64 + if wts != nil { + if wr := len(wts); wr != r { + panic(mat64.ErrShape) + } + + // Weight the rows. + for i := 0; i < c; i++ { + v := xt.RowView(i) + floats.Mul(v, wts) + } + + // Calculate the normalization factor. + n = floats.Sum(wts) + } else { + n = float64(r) + } + + cov.Mul(&xt, &xc) + + // Scale by the sample size. + scale = 1 / (n - 1) + cov.Scale(scale, cov) + + return cov +} diff --git a/covariancematrix_test.go b/covariancematrix_test.go new file mode 100644 index 00000000..63cef50a --- /dev/null +++ b/covariancematrix_test.go @@ -0,0 +1,194 @@ +// Copyright ©2014 The gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package stat + +import ( + "math" + "math/rand" + "testing" + + "github.com/gonum/blas/goblas" + "github.com/gonum/floats" + "github.com/gonum/matrix/mat64" +) + +func init() { + mat64.Register(goblas.Blas{}) +} + +func TestCovarianceMatrix(t *testing.T) { + + // An alternate way to test this is to call the Variance + // and Covariance functions and ensure that the results are identical. + for i, test := range []struct { + data *mat64.Dense + weights mat64.Vec + ans *mat64.Dense + }{ + { + data: mat64.NewDense(5, 2, []float64{ + -2, -4, + -1, 2, + 0, 0, + 1, -2, + 2, 4, + }), + weights: nil, + ans: mat64.NewDense(2, 2, []float64{ + 2.5, 3, + 3, 10, + }), + }, { + data: mat64.NewDense(3, 2, []float64{ + 1, 1, + 2, 4, + 3, 9, + }), + weights: []float64{ + 1, + 1.5, + 1, + }, + ans: mat64.NewDense(2, 2, []float64{ + .8, 3.2, + 3.2, 13.142857142857146, + }), + }, + } { + // Make a copy of the data to check that it isn't changing. + r := test.data.RawMatrix() + d := make([]float64, len(r.Data)) + copy(d, r.Data) + + w := make([]float64, len(test.weights)) + if test.weights != nil { + copy(w, test.weights) + } + c := CovarianceMatrix(nil, test.data, test.weights) + if !c.Equals(test.ans) { + t.Errorf("%d: expected cov %v, found %v", i, test.ans, c) + } + if !floats.Equal(d, r.Data) { + t.Errorf("%d: data was modified during execution") + } + if !floats.Equal(w, test.weights) { + t.Errorf("%d: weights was modified during execution") + } + + // compare with call to Covariance + _, cols := c.Dims() + for ci := 0; ci < cols; ci++ { + for cj := 0; cj < cols; cj++ { + x := test.data.Col(nil, ci) + y := test.data.Col(nil, cj) + cov := Covariance(x, y, test.weights) + if math.Abs(cov-c.At(ci, cj)) > 1e-14 { + t.Errorf("CovMat does not match at (%v, %v). Want %v, got %v.", ci, cj, cov, c.At(ci, cj)) + } + } + } + + } + if !Panics(func() { CovarianceMatrix(nil, mat64.NewDense(5, 2, nil), mat64.Vec([]float64{})) }) { + t.Errorf("CovarianceMatrix did not panic with weight size mismatch") + } + if !Panics(func() { CovarianceMatrix(mat64.NewDense(1, 1, nil), mat64.NewDense(5, 2, nil), nil) }) { + t.Errorf("CovarianceMatrix did not panic with preallocation size mismatch") + } + +} + +// benchmarks + +func randMat(r, c int) mat64.Matrix { + x := make([]float64, r*c) + for i := range x { + x[i] = rand.Float64() + } + return mat64.NewDense(r, c, x) +} + +func benchmarkCovarianceMatrix(b *testing.B, m mat64.Matrix) { + b.ResetTimer() + for i := 0; i < b.N; i++ { + CovarianceMatrix(nil, m, nil) + } +} +func benchmarkCovarianceMatrixInPlace(b *testing.B, m mat64.Matrix) { + _, c := m.Dims() + res := mat64.NewDense(c, c, nil) + b.ResetTimer() + for i := 0; i < b.N; i++ { + CovarianceMatrix(res, m, nil) + } +} + +func BenchmarkCovarianceMatrixSmallxSmall(b *testing.B) { + // 10 * 10 elements + x := randMat(small, small) + benchmarkCovarianceMatrix(b, x) +} +func BenchmarkCovarianceMatrixSmallxMedium(b *testing.B) { + // 10 * 1000 elements + x := randMat(small, medium) + benchmarkCovarianceMatrix(b, x) +} + +func BenchmarkCovarianceMatrixMediumxSmall(b *testing.B) { + // 1000 * 10 elements + x := randMat(medium, small) + benchmarkCovarianceMatrix(b, x) +} +func BenchmarkCovarianceMatrixMediumxMedium(b *testing.B) { + // 1000 * 1000 elements + x := randMat(medium, medium) + benchmarkCovarianceMatrix(b, x) +} + +func BenchmarkCovarianceMatrixLargexSmall(b *testing.B) { + // 1e5 * 10 elements + x := randMat(large, small) + benchmarkCovarianceMatrix(b, x) +} + +func BenchmarkCovarianceMatrixHugexSmall(b *testing.B) { + // 1e7 * 10 elements + x := randMat(huge, small) + benchmarkCovarianceMatrix(b, x) +} + +func BenchmarkCovarianceMatrixSmallxSmallInPlace(b *testing.B) { + // 10 * 10 elements + x := randMat(small, small) + benchmarkCovarianceMatrixInPlace(b, x) +} +func BenchmarkCovarianceMatrixSmallxMediumInPlace(b *testing.B) { + // 10 * 1000 elements + x := randMat(small, medium) + benchmarkCovarianceMatrixInPlace(b, x) +} + +func BenchmarkCovarianceMatrixMediumxSmallInPlace(b *testing.B) { + // 1000 * 10 elements + x := randMat(medium, small) + benchmarkCovarianceMatrixInPlace(b, x) +} +func BenchmarkCovarianceMatrixMediumxMediumInPlace(b *testing.B) { + // 1000 * 1000 elements + x := randMat(medium, medium) + benchmarkCovarianceMatrixInPlace(b, x) +} + +func BenchmarkCovarianceMatrixLargexSmallInPlace(b *testing.B) { + // 1e5 * 10 elements + x := randMat(large, small) + benchmarkCovarianceMatrixInPlace(b, x) +} + +func BenchmarkCovarianceMatrixHugexSmallInPlace(b *testing.B) { + // 1e7 * 10 elements + x := randMat(huge, small) + benchmarkCovarianceMatrixInPlace(b, x) +} diff --git a/moments_bench_test.go b/moments_bench_test.go index 980bf352..6df41394 100644 --- a/moments_bench_test.go +++ b/moments_bench_test.go @@ -17,10 +17,10 @@ import ( ) const ( - SMALL = 10 - MEDIUM = 1000 - LARGE = 100000 - HUGE = 10000000 + small = 10 + medium = 1000 + large = 100000 + huge = 10000000 ) // tests for unweighted versions @@ -41,46 +41,46 @@ func benchmarkMean(b *testing.B, s, wts []float64) { } func BenchmarkMeanSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkMean(b, s, nil) } func BenchmarkMeanMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkMean(b, s, nil) } func BenchmarkMeanLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkMean(b, s, nil) } func BenchmarkMeanHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkMean(b, s, nil) } func BenchmarkMeanSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkMean(b, s, wts) } func BenchmarkMeanMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkMean(b, s, wts) } func BenchmarkMeanLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkMean(b, s, wts) } func BenchmarkMeanHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkMean(b, s, wts) } @@ -92,46 +92,46 @@ func benchmarkVariance(b *testing.B, s, wts []float64) { } func BenchmarkVarianceSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkVariance(b, s, nil) } func BenchmarkVarianceMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkVariance(b, s, nil) } func BenchmarkVarianceLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkVariance(b, s, nil) } func BenchmarkVarianceHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkVariance(b, s, nil) } func BenchmarkVarianceSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkVariance(b, s, wts) } func BenchmarkVarianceMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkVariance(b, s, wts) } func BenchmarkVarianceLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkVariance(b, s, wts) } func BenchmarkVarianceHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkVariance(b, s, wts) } @@ -143,46 +143,46 @@ func benchmarkStdDev(b *testing.B, s, wts []float64) { } func BenchmarkStdDevSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkStdDev(b, s, nil) } func BenchmarkStdDevMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkStdDev(b, s, nil) } func BenchmarkStdDevLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkStdDev(b, s, nil) } func BenchmarkStdDevHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkStdDev(b, s, nil) } func BenchmarkStdDevSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkStdDev(b, s, wts) } func BenchmarkStdDevMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkStdDev(b, s, wts) } func BenchmarkStdDevLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkStdDev(b, s, wts) } func BenchmarkStdDevHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkStdDev(b, s, wts) } @@ -194,46 +194,46 @@ func benchmarkMeanVariance(b *testing.B, s, wts []float64) { } func BenchmarkMeanVarianceSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkMeanVariance(b, s, nil) } func BenchmarkMeanVarianceMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkMeanVariance(b, s, nil) } func BenchmarkMeanVarianceLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkMeanVariance(b, s, nil) } func BenchmarkMeanVarianceHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkMeanVariance(b, s, nil) } func BenchmarkMeanVarianceSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkMeanVariance(b, s, wts) } func BenchmarkMeanVarianceMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkMeanVariance(b, s, wts) } func BenchmarkMeanVarianceLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkMeanVariance(b, s, wts) } func BenchmarkMeanVarianceHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkMeanVariance(b, s, wts) } @@ -245,46 +245,46 @@ func benchmarkMeanStdDev(b *testing.B, s, wts []float64) { } func BenchmarkMeanStdDevSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkMeanStdDev(b, s, nil) } func BenchmarkMeanStdDevMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkMeanStdDev(b, s, nil) } func BenchmarkMeanStdDevLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkMeanStdDev(b, s, nil) } func BenchmarkMeanStdDevHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkMeanStdDev(b, s, nil) } func BenchmarkMeanStdDevSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkMeanStdDev(b, s, wts) } func BenchmarkMeanStdDevMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkMeanStdDev(b, s, wts) } func BenchmarkMeanStdDevLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkMeanStdDev(b, s, wts) } func BenchmarkMeanStdDevHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkMeanStdDev(b, s, wts) } @@ -296,54 +296,54 @@ func benchmarkCovariance(b *testing.B, s1, s2, wts []float64) { } func BenchmarkCovarianceSmall(b *testing.B) { - s1 := RandomSlice(SMALL) - s2 := RandomSlice(SMALL) + s1 := RandomSlice(small) + s2 := RandomSlice(small) benchmarkCovariance(b, s1, s2, nil) } func BenchmarkCovarianceMedium(b *testing.B) { - s1 := RandomSlice(MEDIUM) - s2 := RandomSlice(MEDIUM) + s1 := RandomSlice(medium) + s2 := RandomSlice(medium) benchmarkCovariance(b, s1, s2, nil) } func BenchmarkCovarianceLarge(b *testing.B) { - s1 := RandomSlice(LARGE) - s2 := RandomSlice(LARGE) + s1 := RandomSlice(large) + s2 := RandomSlice(large) benchmarkCovariance(b, s1, s2, nil) } func BenchmarkCovarianceHuge(b *testing.B) { - s1 := RandomSlice(HUGE) - s2 := RandomSlice(HUGE) + s1 := RandomSlice(huge) + s2 := RandomSlice(huge) benchmarkCovariance(b, s1, s2, nil) } func BenchmarkCovarianceSmallWeighted(b *testing.B) { - s1 := RandomSlice(SMALL) - s2 := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s1 := RandomSlice(small) + s2 := RandomSlice(small) + wts := RandomSlice(small) benchmarkCovariance(b, s1, s2, wts) } func BenchmarkCovarianceMediumWeighted(b *testing.B) { - s1 := RandomSlice(MEDIUM) - s2 := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s1 := RandomSlice(medium) + s2 := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkCovariance(b, s1, s2, wts) } func BenchmarkCovarianceLargeWeighted(b *testing.B) { - s1 := RandomSlice(LARGE) - s2 := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s1 := RandomSlice(large) + s2 := RandomSlice(large) + wts := RandomSlice(large) benchmarkCovariance(b, s1, s2, wts) } func BenchmarkCovarianceHugeWeighted(b *testing.B) { - s1 := RandomSlice(HUGE) - s2 := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s1 := RandomSlice(huge) + s2 := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkCovariance(b, s1, s2, wts) } @@ -355,54 +355,54 @@ func benchmarkCorrelation(b *testing.B, s1, s2, wts []float64) { } func BenchmarkCorrelationSmall(b *testing.B) { - s1 := RandomSlice(SMALL) - s2 := RandomSlice(SMALL) + s1 := RandomSlice(small) + s2 := RandomSlice(small) benchmarkCorrelation(b, s1, s2, nil) } func BenchmarkCorrelationMedium(b *testing.B) { - s1 := RandomSlice(MEDIUM) - s2 := RandomSlice(MEDIUM) + s1 := RandomSlice(medium) + s2 := RandomSlice(medium) benchmarkCorrelation(b, s1, s2, nil) } func BenchmarkCorrelationLarge(b *testing.B) { - s1 := RandomSlice(LARGE) - s2 := RandomSlice(LARGE) + s1 := RandomSlice(large) + s2 := RandomSlice(large) benchmarkCorrelation(b, s1, s2, nil) } func BenchmarkCorrelationHuge(b *testing.B) { - s1 := RandomSlice(HUGE) - s2 := RandomSlice(HUGE) + s1 := RandomSlice(huge) + s2 := RandomSlice(huge) benchmarkCorrelation(b, s1, s2, nil) } func BenchmarkCorrelationSmallWeighted(b *testing.B) { - s1 := RandomSlice(SMALL) - s2 := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s1 := RandomSlice(small) + s2 := RandomSlice(small) + wts := RandomSlice(small) benchmarkCorrelation(b, s1, s2, wts) } func BenchmarkCorrelationMediumWeighted(b *testing.B) { - s1 := RandomSlice(MEDIUM) - s2 := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s1 := RandomSlice(medium) + s2 := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkCorrelation(b, s1, s2, wts) } func BenchmarkCorrelationLargeWeighted(b *testing.B) { - s1 := RandomSlice(LARGE) - s2 := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s1 := RandomSlice(large) + s2 := RandomSlice(large) + wts := RandomSlice(large) benchmarkCorrelation(b, s1, s2, wts) } func BenchmarkCorrelationHugeWeighted(b *testing.B) { - s1 := RandomSlice(HUGE) - s2 := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s1 := RandomSlice(huge) + s2 := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkCorrelation(b, s1, s2, wts) } @@ -414,46 +414,46 @@ func benchmarkSkew(b *testing.B, s, wts []float64) { } func BenchmarkSkewSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkSkew(b, s, nil) } func BenchmarkSkewMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkSkew(b, s, nil) } func BenchmarkSkewLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkSkew(b, s, nil) } func BenchmarkSkewHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkSkew(b, s, nil) } func BenchmarkSkewSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkSkew(b, s, wts) } func BenchmarkSkewMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkSkew(b, s, wts) } func BenchmarkSkewLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkSkew(b, s, wts) } func BenchmarkSkewHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkSkew(b, s, wts) } @@ -465,46 +465,46 @@ func benchmarkExKurtosis(b *testing.B, s, wts []float64) { } func BenchmarkExKurtosisSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkExKurtosis(b, s, nil) } func BenchmarkExKurtosisMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkExKurtosis(b, s, nil) } func BenchmarkExKurtosisLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkExKurtosis(b, s, nil) } func BenchmarkExKurtosisHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkExKurtosis(b, s, nil) } func BenchmarkExKurtosisSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkExKurtosis(b, s, wts) } func BenchmarkExKurtosisMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkExKurtosis(b, s, wts) } func BenchmarkExKurtosisLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkExKurtosis(b, s, wts) } func BenchmarkExKurtosisHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkExKurtosis(b, s, wts) } @@ -516,46 +516,46 @@ func benchmarkMoment(b *testing.B, n float64, s, wts []float64) { } func BenchmarkMomentSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkMoment(b, 5, s, nil) } func BenchmarkMomentMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkMoment(b, 5, s, nil) } func BenchmarkMomentLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkMoment(b, 5, s, nil) } func BenchmarkMomentHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkMoment(b, 5, s, nil) } func BenchmarkMomentSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkMoment(b, 5, s, wts) } func BenchmarkMomentMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkMoment(b, 5, s, wts) } func BenchmarkMomentLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkMoment(b, 5, s, wts) } func BenchmarkMomentHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkMoment(b, 5, s, wts) } @@ -567,45 +567,45 @@ func benchmarkMomentAbout(b *testing.B, n float64, s []float64, mean float64, wt } func BenchmarkMomentAboutSmall(b *testing.B) { - s := RandomSlice(SMALL) + s := RandomSlice(small) benchmarkMomentAbout(b, 5, s, 0, nil) } func BenchmarkMomentAboutMedium(b *testing.B) { - s := RandomSlice(MEDIUM) + s := RandomSlice(medium) benchmarkMomentAbout(b, 5, s, 0, nil) } func BenchmarkMomentAboutLarge(b *testing.B) { - s := RandomSlice(LARGE) + s := RandomSlice(large) benchmarkMomentAbout(b, 5, s, 0, nil) } func BenchmarkMomentAboutHuge(b *testing.B) { - s := RandomSlice(HUGE) + s := RandomSlice(huge) benchmarkMomentAbout(b, 5, s, 0, nil) } func BenchmarkMomentAboutSmallWeighted(b *testing.B) { - s := RandomSlice(SMALL) - wts := RandomSlice(SMALL) + s := RandomSlice(small) + wts := RandomSlice(small) benchmarkMomentAbout(b, 5, s, 0, wts) } func BenchmarkMomentAboutMediumWeighted(b *testing.B) { - s := RandomSlice(MEDIUM) - wts := RandomSlice(MEDIUM) + s := RandomSlice(medium) + wts := RandomSlice(medium) benchmarkMomentAbout(b, 5, s, 0, wts) } func BenchmarkMomentAboutLargeWeighted(b *testing.B) { - s := RandomSlice(LARGE) - wts := RandomSlice(LARGE) + s := RandomSlice(large) + wts := RandomSlice(large) benchmarkMomentAbout(b, 5, s, 0, wts) } func BenchmarkMomentAboutHugeWeighted(b *testing.B) { - s := RandomSlice(HUGE) - wts := RandomSlice(HUGE) + s := RandomSlice(huge) + wts := RandomSlice(huge) benchmarkMomentAbout(b, 5, s, 0, wts) } diff --git a/stat.go b/stat.go index 1d616af1..4c76fbb3 100644 --- a/stat.go +++ b/stat.go @@ -232,10 +232,10 @@ func Covariance(x, y, weights []float64) float64 { w := weights[i] yv := y[i] wxd := w * (xv - xu) - wyd := w * (yv - yu) - ss += wxd * wyd + yd := (yv - yu) + ss += wxd * yd xcompensation += wxd - ycompensation += wyd + ycompensation += w * yd sumWeights += w } // xcompensation and ycompensation are from Chan, et. al. diff --git a/stat_test.go b/stat_test.go index c646aad9..51e9a54c 100644 --- a/stat_test.go +++ b/stat_test.go @@ -128,6 +128,12 @@ func TestCovariance(t *testing.T) { weights: []float64{1, 1.5, 1}, ans: 3.2, }, + { + p: []float64{1, 4, 9}, + q: []float64{1, 4, 9}, + weights: []float64{1, 1.5, 1}, + ans: 13.142857142857146, + }, } { c := Covariance(test.p, test.q, test.weights) if math.Abs(c-test.ans) > 1e-14 { @@ -1256,6 +1262,16 @@ func TestVariance(t *testing.T) { weights: []float64{2, 1, 2, 1, 1}, ans: 4.2857142857142865, }, + { + x: []float64{1, 4, 9}, + weights: []float64{1, 1.5, 1}, + ans: 13.142857142857146, + }, + { + x: []float64{1, 2, 3}, + weights: []float64{1, 1.5, 1}, + ans: .8, + }, } { variance := Variance(test.x, test.weights) if math.Abs(variance-test.ans) > 1e-14 {