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
|
// Copyright ©2016 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 mathext
import (
"testing"
"gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/floats/scalar"
)
func TestIncBeta(t *testing.T) {
t.Parallel()
tol := 1e-14
tol2 := 1e-10
// Test against values from scipy
for i, test := range []struct {
a, b, x, ans float64
}{
{1, 1, 0.8, 0.8},
{1, 5, 0.8, 0.99968000000000001},
{10, 10, 0.8, 0.99842087945083291},
{10, 10, 0.1, 3.929882327128003e-06},
{10, 2, 0.4, 0.00073400320000000028},
{0.1, 0.2, 0.6, 0.69285678232066683},
{1, 10, 0.7489, 0.99999900352334858},
} {
y := RegIncBeta(test.a, test.b, test.x)
if !scalar.EqualWithinAbsOrRel(y, test.ans, tol, tol) {
t.Errorf("Incomplete beta mismatch. Case %v: Got %v, want %v", i, y, test.ans)
}
yc := 1 - RegIncBeta(test.b, test.a, 1-test.x)
if !scalar.EqualWithinAbsOrRel(y, yc, tol, tol) {
t.Errorf("Incomplete beta complementary mismatch. Case %v: Got %v, want %v", i, y, yc)
}
x := InvRegIncBeta(test.a, test.b, y)
if !scalar.EqualWithinAbsOrRel(x, test.x, tol2, tol2) {
t.Errorf("Inverse incomplete beta mismatch. Case %v: Got %v, want %v", i, x, test.x)
}
}
// Confirm that Invincbeta and Incbeta agree. Sweep over a variety of
// a, b, and y values.
tol = 1e-6
steps := 201
ints := make([]float64, steps)
floats.Span(ints, 0, 1)
sz := 51
min := 1e-2
max := 1e2
as := make([]float64, sz)
floats.LogSpan(as, min, max)
bs := make([]float64, sz)
floats.LogSpan(bs, min, max)
for _, a := range as {
for _, b := range bs {
for _, yr := range ints {
x := InvRegIncBeta(a, b, yr)
if x > 1-1e-6 {
// Numerical error too large
continue
}
y := RegIncBeta(a, b, x)
if !scalar.EqualWithinAbsOrRel(yr, y, tol, tol) {
t.Errorf("Mismatch between inv inc beta and inc beta. a = %v, b = %v, x = %v, got %v, want %v.", a, b, x, y, yr)
break
}
}
}
}
}
|
