aboutsummaryrefslogtreecommitdiff
path: root/vendor/github.com/aclements/go-moremath/mathx/beta.go
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/github.com/aclements/go-moremath/mathx/beta.go')
-rw-r--r--vendor/github.com/aclements/go-moremath/mathx/beta.go93
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")
-}