github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/mathext/beta_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_test
     6  
     7  import (
     8  	"math"
     9  	"testing"
    10  
    11  	"github.com/jingcheng-WU/gonum/floats/scalar"
    12  	"github.com/jingcheng-WU/gonum/mathext"
    13  )
    14  
    15  var betaTests = []struct {
    16  	p, q float64
    17  	want float64
    18  }{
    19  	{
    20  		p:    1,
    21  		q:    2,
    22  		want: 0.5, // obtained from scipy.special.beta(1,2) (version=0.18.0)
    23  	},
    24  	{
    25  		p:    10,
    26  		q:    20,
    27  		want: 4.9925087406346778e-09, // obtained from scipy.special.beta(10,20) (version=0.18.0)
    28  	},
    29  	{
    30  		p:    +0,
    31  		q:    10,
    32  		want: math.Inf(+1),
    33  	},
    34  	{
    35  		p:    -0,
    36  		q:    10,
    37  		want: math.Inf(+1),
    38  	},
    39  	{
    40  		p:    0,
    41  		q:    0,
    42  		want: math.NaN(),
    43  	},
    44  	{
    45  		p:    0,
    46  		q:    math.Inf(-1),
    47  		want: math.NaN(),
    48  	},
    49  	{
    50  		p:    10,
    51  		q:    math.Inf(-1),
    52  		want: math.NaN(),
    53  	},
    54  	{
    55  		p:    0,
    56  		q:    math.Inf(+1),
    57  		want: math.NaN(),
    58  	},
    59  	{
    60  		p:    10,
    61  		q:    math.Inf(+1),
    62  		want: math.NaN(),
    63  	},
    64  	{
    65  		p:    math.NaN(),
    66  		q:    10,
    67  		want: math.NaN(),
    68  	},
    69  	{
    70  		p:    math.NaN(),
    71  		q:    0,
    72  		want: math.NaN(),
    73  	},
    74  	{
    75  		p:    -1,
    76  		q:    0,
    77  		want: math.NaN(),
    78  	},
    79  	{
    80  		p:    -1,
    81  		q:    +1,
    82  		want: math.NaN(),
    83  	},
    84  }
    85  
    86  func TestBeta(t *testing.T) {
    87  	t.Parallel()
    88  	for i, test := range betaTests {
    89  		v := mathext.Beta(test.p, test.q)
    90  		testOK := func(x float64) bool {
    91  			return scalar.EqualWithinAbsOrRel(x, test.want, 1e-15, 1e-15) || (math.IsNaN(test.want) && math.IsNaN(x))
    92  		}
    93  		if !testOK(v) {
    94  			t.Errorf("test #%d: Beta(%v, %v)=%v. want=%v\n",
    95  				i, test.p, test.q, v, test.want,
    96  			)
    97  		}
    98  
    99  		u := mathext.Beta(test.q, test.p)
   100  		if !testOK(u) {
   101  			t.Errorf("test #%[1]d: Beta(%[2]v, %[3]v)=%[4]v != Beta(%[3]v, %[2]v)=%[5]v)\n",
   102  				i, test.p, test.q, v, u,
   103  			)
   104  		}
   105  
   106  		if math.IsInf(v, +1) || math.IsNaN(v) {
   107  			continue
   108  		}
   109  
   110  		vv := mathext.Beta(test.p, test.q+1)
   111  		uu := mathext.Beta(test.p+1, test.q)
   112  		if !scalar.EqualWithinAbsOrRel(v, vv+uu, 1e-15, 1e-15) {
   113  			t.Errorf(
   114  				"test #%[1]d: Beta(%[2]v, %[3]v)=%[4]v != Beta(%[2]v+1, %[3]v) + Beta(%[2]v, %[3]v+1) (=%[5]v + %[6]v = %[7]v)\n",
   115  				i, test.p, test.q, v, uu, vv, uu+vv,
   116  			)
   117  		}
   118  
   119  		vbeta2 := beta2(test.p, test.q)
   120  		if !scalar.EqualWithinAbsOrRel(v, vbeta2, 1e-15, 1e-15) {
   121  			t.Errorf(
   122  				"test #%[1]d: Beta(%[2]v, %[3]v) != Γ(p)Γ(q) / Γ(p+q) (v=%[4]v u=%[5]v)\n",
   123  				i, test.p, test.q, v, vbeta2,
   124  			)
   125  		}
   126  	}
   127  }
   128  
   129  func beta2(x, y float64) float64 {
   130  	return math.Gamma(x) * math.Gamma(y) / math.Gamma(x+y)
   131  }
   132  
   133  func BenchmarkBeta(b *testing.B) {
   134  	for i := 0; i < b.N; i++ {
   135  		_ = mathext.Beta(10, 20)
   136  	}
   137  }
   138  
   139  func BenchmarkBeta2(b *testing.B) {
   140  	for i := 0; i < b.N; i++ {
   141  		_ = math.Gamma(10) * math.Gamma(20) / math.Gamma(10+20)
   142  	}
   143  }
   144  
   145  func TestLbeta(t *testing.T) {
   146  	t.Parallel()
   147  	for i, test := range betaTests {
   148  		want := math.Log(test.want)
   149  		v := mathext.Lbeta(test.p, test.q)
   150  
   151  		testOK := func(x float64) bool {
   152  			return scalar.EqualWithinAbsOrRel(x, want, 1e-15, 1e-15) || (math.IsNaN(want) && math.IsNaN(x))
   153  		}
   154  		if !testOK(v) {
   155  			t.Errorf("test #%d: Lbeta(%v, %v)=%v. want=%v\n",
   156  				i, test.p, test.q, v, want,
   157  			)
   158  		}
   159  
   160  		u := mathext.Lbeta(test.q, test.p)
   161  		if !testOK(u) {
   162  			t.Errorf("test #%[1]d: Lbeta(%[2]v, %[3]v)=%[4]v != Lbeta(%[3]v, %[2]v)=%[5]v)\n",
   163  				i, test.p, test.q, v, u,
   164  			)
   165  		}
   166  
   167  		if math.IsInf(v, +1) || math.IsNaN(v) {
   168  			continue
   169  		}
   170  
   171  		vbeta2 := math.Log(beta2(test.p, test.q))
   172  		if !scalar.EqualWithinAbsOrRel(v, vbeta2, 1e-15, 1e-15) {
   173  			t.Errorf(
   174  				"test #%[1]d: Lbeta(%[2]v, %[3]v) != Log(Γ(p)Γ(q) / Γ(p+q)) (v=%[4]v u=%[5]v)\n",
   175  				i, test.p, test.q, v, vbeta2,
   176  			)
   177  		}
   178  	}
   179  }