diff options
Diffstat (limited to 'vendor/github.com/aclements/go-moremath/mathx/beta.go')
| -rw-r--r-- | vendor/github.com/aclements/go-moremath/mathx/beta.go | 93 |
1 files changed, 0 insertions, 93 deletions
diff --git a/vendor/github.com/aclements/go-moremath/mathx/beta.go b/vendor/github.com/aclements/go-moremath/mathx/beta.go deleted file mode 100644 index 49f8722..0000000 --- a/vendor/github.com/aclements/go-moremath/mathx/beta.go +++ /dev/null @@ -1,93 +0,0 @@ -// Copyright 2015 The Go 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 mathx - -import "math" - -func lgamma(x float64) float64 { - y, _ := math.Lgamma(x) - return y -} - -// Beta returns the value of the complete beta function B(a, b). -func Beta(a, b float64) float64 { - // B(x,y) = Γ(x)Γ(y) / Γ(x+y) - return math.Exp(lgamma(a) + lgamma(b) - lgamma(a+b)) -} - -// BetaInc returns the value of the regularized incomplete beta -// function Iₓ(a, b). -// -// This is not to be confused with the "incomplete beta function", -// which can be computed as BetaInc(x, a, b)*Beta(a, b). -// -// If x < 0 or x > 1, returns NaN. -func BetaInc(x, a, b float64) float64 { - // Based on Numerical Recipes in C, section 6.4. This uses the - // continued fraction definition of I: - // - // (xᵃ*(1-x)ᵇ)/(a*B(a,b)) * (1/(1+(d₁/(1+(d₂/(1+...)))))) - // - // where B(a,b) is the beta function and - // - // d_{2m+1} = -(a+m)(a+b+m)x/((a+2m)(a+2m+1)) - // d_{2m} = m(b-m)x/((a+2m-1)(a+2m)) - if x < 0 || x > 1 { - return math.NaN() - } - bt := 0.0 - if 0 < x && x < 1 { - // Compute the coefficient before the continued - // fraction. - bt = math.Exp(lgamma(a+b) - lgamma(a) - lgamma(b) + - a*math.Log(x) + b*math.Log(1-x)) - } - if x < (a+1)/(a+b+2) { - // Compute continued fraction directly. - return bt * betacf(x, a, b) / a - } else { - // Compute continued fraction after symmetry transform. - return 1 - bt*betacf(1-x, b, a)/b - } -} - -// betacf is the continued fraction component of the regularized -// incomplete beta function Iₓ(a, b). -func betacf(x, a, b float64) float64 { - const maxIterations = 200 - const epsilon = 3e-14 - - raiseZero := func(z float64) float64 { - if math.Abs(z) < math.SmallestNonzeroFloat64 { - return math.SmallestNonzeroFloat64 - } - return z - } - - c := 1.0 - d := 1 / raiseZero(1-(a+b)*x/(a+1)) - h := d - for m := 1; m <= maxIterations; m++ { - mf := float64(m) - - // Even step of the recurrence. - numer := mf * (b - mf) * x / ((a + 2*mf - 1) * (a + 2*mf)) - d = 1 / raiseZero(1+numer*d) - c = raiseZero(1 + numer/c) - h *= d * c - - // Odd step of the recurrence. - numer = -(a + mf) * (a + b + mf) * x / ((a + 2*mf) * (a + 2*mf + 1)) - d = 1 / raiseZero(1+numer*d) - c = raiseZero(1 + numer/c) - hfac := d * c - h *= hfac - - if math.Abs(hfac-1) < epsilon { - return h - } - } - panic("betainc: a or b too big; failed to converge") -} |
