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  }