Skip to content

Commit 40090d8

Browse files
ncrucesncruces
committed
Moments.
1 parent d2f1629 commit 40090d8

File tree

6 files changed

+202
-79
lines changed

6 files changed

+202
-79
lines changed

‎ext/stats/kahan.go

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
package stats
2+
3+
// https://en.wikipedia.org/wiki/Kahan_summation_algorithm
4+
5+
type kahan struct{ hi, lo float64 }
6+
7+
func (k *kahan) add(x float64) {
8+
y := k.lo + x
9+
t := k.hi + y
10+
k.lo = y - (t - k.hi)
11+
k.hi = t
12+
}
13+
14+
func (k *kahan) sub(x float64) {
15+
y := k.lo - x
16+
t := k.hi + y
17+
k.lo = y - (t - k.hi)
18+
k.hi = t
19+
}

‎ext/stats/moments.go

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
package stats
2+
3+
import "math"
4+
5+
// Terriberry's algorithm with Kahan summation:
6+
// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Higher-order_statistics
7+
8+
type moments struct {
9+
m1, m2, m3, m4 kahan
10+
n int64
11+
}
12+
13+
func (m moments) mean() float64 {
14+
return m.m1.hi
15+
}
16+
17+
func (m moments) var_pop() float64 {
18+
return m.m2.hi / float64(m.n)
19+
}
20+
21+
func (m moments) var_samp() float64 {
22+
return m.m2.hi / float64(m.n-1) // Bessel's correction
23+
}
24+
25+
func (m moments) stddev_pop() float64 {
26+
return math.Sqrt(m.var_pop())
27+
}
28+
29+
func (m moments) stddev_samp() float64 {
30+
return math.Sqrt(m.var_samp())
31+
}
32+
33+
func (m moments) skewness_pop() float64 {
34+
m2 := m.m2.hi
35+
if div := m2 * m2 * m2; div != 0 {
36+
return m.m3.hi * math.Sqrt(float64(m.n)/div)
37+
}
38+
return math.NaN()
39+
}
40+
41+
func (m moments) skewness_samp() float64 {
42+
n := m.n
43+
// https://mathworks.com/help/stats/skewness.html#f1132178
44+
return m.skewness_pop() * math.Sqrt(float64(n*(n-1))) / float64(n-2)
45+
}
46+
47+
func (m moments) kurtosis_pop() float64 {
48+
m2 := m.m2.hi
49+
if div := m2 * m2; div != 0 {
50+
return m.m4.hi * float64(m.n) / div
51+
}
52+
return math.NaN()
53+
}
54+
55+
func (m moments) kurtosis_samp() float64 {
56+
n := m.n
57+
// https://mathworks.com/help/stats/kurtosis.html#f4975293
58+
k := math.FMA(m.kurtosis_pop(), float64(n+1), float64(3-3*n))
59+
return math.FMA(k, float64(n-1)/float64((n-2)*(n-3)), 3)
60+
}
61+
62+
func (m *moments) enqueue(x float64) {
63+
n := m.n + 1
64+
m.n = n
65+
d1 := x - m.m1.hi - m.m1.lo
66+
dn := d1 / float64(n)
67+
d2 := dn * dn
68+
t1 := d1 * dn * float64(n-1)
69+
m.m4.add(t1*d2*float64(n*n-3*n+3) + 6*d2*m.m2.hi - 4*dn*m.m3.hi)
70+
m.m3.add(t1*dn*float64(n-2) - 3*dn*m.m2.hi)
71+
m.m2.add(t1)
72+
m.m1.add(dn)
73+
}
74+
75+
func (m *moments) dequeue(x float64) {
76+
n := m.n - 1
77+
if n <= 0 {
78+
*m = moments{}
79+
return
80+
}
81+
m.n = n
82+
d1 := x - m.m1.hi - m.m1.lo
83+
dn := d1 / float64(n)
84+
d2 := dn * dn
85+
t1 := d1 * dn * float64(n+1)
86+
m.m4.sub(t1*d2*float64(n*n+3*n+3) - 6*d2*m.m2.hi - 4*dn*m.m3.hi)
87+
m.m3.sub(t1*dn*float64(n+2) - 3*dn*m.m2.hi)
88+
m.m2.sub(t1)
89+
m.m1.sub(dn)
90+
}

‎ext/stats/moments.nogo

Lines changed: 0 additions & 38 deletions
This file was deleted.

‎ext/stats/moments_test.go

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
package stats
2+
3+
import (
4+
"math"
5+
"testing"
6+
)
7+
8+
func Test_moments(t *testing.T) {
9+
t.Parallel()
10+
11+
var s1 moments
12+
s1.enqueue(1)
13+
s1.dequeue(1)
14+
if !math.IsNaN(s1.skewness_pop()) {
15+
t.Errorf("want NaN")
16+
}
17+
if !math.IsNaN(s1.kurtosis_pop()) {
18+
t.Errorf("want NaN")
19+
}
20+
21+
s1.enqueue(+0.5377)
22+
s1.enqueue(+1.8339)
23+
s1.enqueue(-2.2588)
24+
s1.enqueue(+0.8622)
25+
s1.enqueue(+0.3188)
26+
s1.enqueue(-1.3077)
27+
s1.enqueue(-0.4336)
28+
s1.enqueue(+0.3426)
29+
s1.enqueue(+3.5784)
30+
s1.enqueue(+2.7694)
31+
32+
if got := float32(s1.skewness_pop()); got != 0.106098293 {
33+
t.Errorf("got %v, want 0.1061", got)
34+
}
35+
if got := float32(s1.skewness_samp()); got != 0.1258171 {
36+
t.Errorf("got %v, want 0.1258", got)
37+
}
38+
if got := float32(s1.kurtosis_pop()); got != 2.3121266 {
39+
t.Errorf("got %v, want 2.3121", got)
40+
}
41+
if got := float32(s1.kurtosis_samp()); got != 2.7482237 {
42+
t.Errorf("got %v, want 2.7483", got)
43+
}
44+
45+
var s2 welford
46+
47+
s2.enqueue(+0.5377)
48+
s2.enqueue(+1.8339)
49+
s2.enqueue(-2.2588)
50+
s2.enqueue(+0.8622)
51+
s2.enqueue(+0.3188)
52+
s2.enqueue(-1.3077)
53+
s2.enqueue(-0.4336)
54+
s2.enqueue(+0.3426)
55+
s2.enqueue(+3.5784)
56+
s2.enqueue(+2.7694)
57+
58+
if got, want := s1.mean(), s2.mean(); got != want {
59+
t.Errorf("got %v, want %v", got, want)
60+
}
61+
if got, want := s1.stddev_pop(), s2.stddev_pop(); got != want {
62+
t.Errorf("got %v, want %v", got, want)
63+
}
64+
if got, want := s1.stddev_samp(), s2.stddev_samp(); got != want {
65+
t.Errorf("got %v, want %v", got, want)
66+
}
67+
68+
s1.enqueue(math.Pi)
69+
s1.enqueue(math.Sqrt2)
70+
s1.enqueue(math.E)
71+
s1.dequeue(math.Pi)
72+
s1.dequeue(math.E)
73+
s1.dequeue(math.Sqrt2)
74+
75+
if got := float32(s1.skewness_pop()); got != 0.106098293 {
76+
t.Errorf("got %v, want 0.1061", got)
77+
}
78+
if got := float32(s1.skewness_samp()); got != 0.1258171 {
79+
t.Errorf("got %v, want 0.1258", got)
80+
}
81+
if got := float32(s1.kurtosis_pop()); got != 2.3121266 {
82+
t.Errorf("got %v, want 2.3121", got)
83+
}
84+
if got := float32(s1.kurtosis_samp()); got != 2.7482237 {
85+
t.Errorf("got %v, want 2.7483", got)
86+
}
87+
}

‎ext/stats/welford.go

Lines changed: 1 addition & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -10,14 +10,13 @@ import (
1010
// Welford's algorithm with Kahan summation:
1111
// The effect of truncation in statistical computation [van Reeken, AJ 1970]
1212
// https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm
13-
// https://en.wikipedia.org/wiki/Kahan_summation_algorithm
1413

1514
type welford struct {
1615
m1, m2 kahan
1716
n int64
1817
}
1918

20-
func (w welford) average() float64 {
19+
func (w welford) mean() float64 {
2120
return w.m1.hi
2221
}
2322

@@ -170,19 +169,3 @@ func (w *welford2) dequeue(y, x float64) {
170169
w.m2x.sub(d1x * d2x)
171170
w.cov.sub(d1y * d2x)
172171
}
173-
174-
type kahan struct{ hi, lo float64 }
175-
176-
func (k *kahan) add(x float64) {
177-
y := k.lo + x
178-
t := k.hi + y
179-
k.lo = y - (t - k.hi)
180-
k.hi = t
181-
}
182-
183-
func (k *kahan) sub(x float64) {
184-
y := k.lo - x
185-
t := k.hi + y
186-
k.lo = y - (t - k.hi)
187-
k.hi = t
188-
}

‎ext/stats/welford_test.go

Lines changed: 5 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -9,12 +9,14 @@ func Test_welford(t *testing.T) {
99
t.Parallel()
1010

1111
var s1, s2 welford
12+
s1.enqueue(1)
13+
s1.dequeue(1)
1214

1315
s1.enqueue(4)
1416
s1.enqueue(7)
1517
s1.enqueue(13)
1618
s1.enqueue(16)
17-
if got := s1.average(); got != 10 {
19+
if got := s1.mean(); got != 10 {
1820
t.Errorf("got %v, want 10", got)
1921
}
2022
if got := s1.var_samp(); got != 30 {
@@ -37,22 +39,14 @@ func Test_welford(t *testing.T) {
3739
if s1.var_pop() != s2.var_pop() {
3840
t.Errorf("got %v, want %v", s1, s2)
3941
}
40-
41-
s1.dequeue(16)
42-
s1.dequeue(7)
43-
s1.dequeue(13)
44-
s1.enqueue(16)
45-
s1.enqueue(7)
46-
s1.enqueue(13)
47-
if s1.var_pop() != s2.var_pop() {
48-
t.Errorf("got %v, want %v", s1, s2)
49-
}
5042
}
5143

5244
func Test_covar(t *testing.T) {
5345
t.Parallel()
5446

5547
var c1, c2 welford2
48+
c1.enqueue(1, 1)
49+
c1.dequeue(1, 1)
5650

5751
c1.enqueue(3, 70)
5852
c1.enqueue(5, 80)
@@ -75,18 +69,6 @@ func Test_covar(t *testing.T) {
7569
if c1.covar_pop() != c2.covar_pop() {
7670
t.Errorf("got %v, want %v", c1.covar_pop(), c2.covar_pop())
7771
}
78-
79-
c1.dequeue(2, 60)
80-
c1.dequeue(5, 80)
81-
c1.dequeue(4, 75)
82-
c1.dequeue(7, 90)
83-
c1.enqueue(2, 60)
84-
c1.enqueue(5, 80)
85-
c1.enqueue(4, 75)
86-
c1.enqueue(7, 90)
87-
if c1.covar_pop() != c2.covar_pop() {
88-
t.Errorf("got %v, want %v", c1.covar_pop(), c2.covar_pop())
89-
}
9072
}
9173

9274
func Test_correlation(t *testing.T) {

0 commit comments

Comments
 (0)