github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/betainc_test.go (about)

     1  // Copyright ©2016 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package mathext
     6  
     7  import (
     8  	"testing"
     9  
    10  	"github.com/jingcheng-WU/gonum/floats"
    11  	"github.com/jingcheng-WU/gonum/floats/scalar"
    12  )
    13  
    14  func TestIncBeta(t *testing.T) {
    15  	t.Parallel()
    16  	tol := 1e-14
    17  	tol2 := 1e-10
    18  	// Test against values from scipy
    19  	for i, test := range []struct {
    20  		a, b, x, ans float64
    21  	}{
    22  		{1, 1, 0.8, 0.8},
    23  		{1, 5, 0.8, 0.99968000000000001},
    24  		{10, 10, 0.8, 0.99842087945083291},
    25  		{10, 10, 0.1, 3.929882327128003e-06},
    26  		{10, 2, 0.4, 0.00073400320000000028},
    27  		{0.1, 0.2, 0.6, 0.69285678232066683},
    28  		{1, 10, 0.7489, 0.99999900352334858},
    29  	} {
    30  		y := RegIncBeta(test.a, test.b, test.x)
    31  		if !scalar.EqualWithinAbsOrRel(y, test.ans, tol, tol) {
    32  			t.Errorf("Incomplete beta mismatch. Case %v: Got %v, want %v", i, y, test.ans)
    33  		}
    34  
    35  		yc := 1 - RegIncBeta(test.b, test.a, 1-test.x)
    36  		if !scalar.EqualWithinAbsOrRel(y, yc, tol, tol) {
    37  			t.Errorf("Incomplete beta complementary mismatch. Case %v: Got %v, want %v", i, y, yc)
    38  		}
    39  
    40  		x := InvRegIncBeta(test.a, test.b, y)
    41  		if !scalar.EqualWithinAbsOrRel(x, test.x, tol2, tol2) {
    42  			t.Errorf("Inverse incomplete beta mismatch. Case %v: Got %v, want %v", i, x, test.x)
    43  		}
    44  	}
    45  
    46  	// Confirm that Invincbeta and Incbeta agree. Sweep over a variety of
    47  	// a, b, and y values.
    48  	tol = 1e-6
    49  	steps := 201
    50  	ints := make([]float64, steps)
    51  	floats.Span(ints, 0, 1)
    52  
    53  	sz := 51
    54  	min := 1e-2
    55  	max := 1e2
    56  	as := make([]float64, sz)
    57  	floats.LogSpan(as, min, max)
    58  	bs := make([]float64, sz)
    59  	floats.LogSpan(bs, min, max)
    60  
    61  	for _, a := range as {
    62  		for _, b := range bs {
    63  			for _, yr := range ints {
    64  				x := InvRegIncBeta(a, b, yr)
    65  				if x > 1-1e-6 {
    66  					// Numerical error too large
    67  					continue
    68  				}
    69  				y := RegIncBeta(a, b, x)
    70  				if !scalar.EqualWithinAbsOrRel(yr, y, tol, tol) {
    71  					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)
    72  					break
    73  				}
    74  			}
    75  		}
    76  	}
    77  }