github.com/jgbaldwinbrown/perf@v0.1.1/pkg/stats/beta.go (about) 1 // Copyright 2015 The Go 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 stats 6 7 import "math" 8 9 func lgamma(x float64) float64 { 10 y, _ := math.Lgamma(x) 11 return y 12 } 13 14 // mathBeta returns the value of the complete beta function B(a, b). 15 func mathBeta(a, b float64) float64 { 16 // B(x,y) = Γ(x)Γ(y) / Γ(x+y) 17 return math.Exp(lgamma(a) + lgamma(b) - lgamma(a+b)) 18 } 19 20 // mathBetaInc returns the value of the regularized incomplete beta 21 // function Iₓ(a, b). 22 // 23 // This is not to be confused with the "incomplete beta function", 24 // which can be computed as BetaInc(x, a, b)*Beta(a, b). 25 // 26 // If x < 0 or x > 1, returns NaN. 27 func mathBetaInc(x, a, b float64) float64 { 28 // Based on Numerical Recipes in C, section 6.4. This uses the 29 // continued fraction definition of I: 30 // 31 // (xᵃ*(1-x)ᵇ)/(a*B(a,b)) * (1/(1+(d₁/(1+(d₂/(1+...)))))) 32 // 33 // where B(a,b) is the beta function and 34 // 35 // d_{2m+1} = -(a+m)(a+b+m)x/((a+2m)(a+2m+1)) 36 // d_{2m} = m(b-m)x/((a+2m-1)(a+2m)) 37 if x < 0 || x > 1 { 38 return math.NaN() 39 } 40 bt := 0.0 41 if 0 < x && x < 1 { 42 // Compute the coefficient before the continued 43 // fraction. 44 bt = math.Exp(lgamma(a+b) - lgamma(a) - lgamma(b) + 45 a*math.Log(x) + b*math.Log(1-x)) 46 } 47 if x < (a+1)/(a+b+2) { 48 // Compute continued fraction directly. 49 return bt * betacf(x, a, b) / a 50 } else { 51 // Compute continued fraction after symmetry transform. 52 return 1 - bt*betacf(1-x, b, a)/b 53 } 54 } 55 56 // betacf is the continued fraction component of the regularized 57 // incomplete beta function Iₓ(a, b). 58 func betacf(x, a, b float64) float64 { 59 const maxIterations = 200 60 const epsilon = 3e-14 61 62 raiseZero := func(z float64) float64 { 63 if math.Abs(z) < math.SmallestNonzeroFloat64 { 64 return math.SmallestNonzeroFloat64 65 } 66 return z 67 } 68 69 c := 1.0 70 d := 1 / raiseZero(1-(a+b)*x/(a+1)) 71 h := d 72 for m := 1; m <= maxIterations; m++ { 73 mf := float64(m) 74 75 // Even step of the recurrence. 76 numer := mf * (b - mf) * x / ((a + 2*mf - 1) * (a + 2*mf)) 77 d = 1 / raiseZero(1+numer*d) 78 c = raiseZero(1 + numer/c) 79 h *= d * c 80 81 // Odd step of the recurrence. 82 numer = -(a + mf) * (a + b + mf) * x / ((a + 2*mf) * (a + 2*mf + 1)) 83 d = 1 / raiseZero(1+numer*d) 84 c = raiseZero(1 + numer/c) 85 hfac := d * c 86 h *= hfac 87 88 if math.Abs(hfac-1) < epsilon { 89 return h 90 } 91 } 92 panic("betainc: a or b too big; failed to converge") 93 }