gonum.org/v1/gonum@v0.14.0/mathext/gamma_inc_inv.go (about)

     1  // Derived from SciPy's special/c_misc/gammaincinv.c
     2  // https://github.com/scipy/scipy/blob/master/scipy/special/c_misc/gammaincinv.c
     3  
     4  // Copyright ©2017 The Gonum Authors. All rights reserved.
     5  // Use of this source code is governed by a BSD-style
     6  // license that can be found in the LICENSE file.
     7  
     8  package mathext
     9  
    10  import (
    11  	"math"
    12  
    13  	"gonum.org/v1/gonum/mathext/internal/cephes"
    14  )
    15  
    16  const (
    17  	allowedATol = 1e-306
    18  	allowedRTol = 1e-6
    19  )
    20  
    21  func gammaIncReg(x float64, params []float64) float64 {
    22  	return cephes.Igam(params[0], x) - params[1]
    23  }
    24  
    25  // gammaIncRegInv is the inverse of the regularized incomplete Gamma integral. That is, it
    26  // returns x such that:
    27  //
    28  //	Igam(a, x) = y
    29  //
    30  // The input argument a must be positive and y must be between 0 and 1
    31  // inclusive or gammaIncRegInv will panic. gammaIncRegInv should return a
    32  // positive number, but can return NaN if there is a failure to converge.
    33  func gammaIncRegInv(a, y float64) float64 {
    34  	// For y not small, we just use
    35  	//  IgamI(a, 1-y)
    36  	// (inverse of the complemented incomplete Gamma integral). For y small,
    37  	// however, 1-y is about 1, and we lose digits.
    38  	if a <= 0 || y <= 0 || y >= 0.25 {
    39  		return cephes.IgamI(a, 1-y)
    40  	}
    41  
    42  	lo := 0.0
    43  	flo := -y
    44  	hi := cephes.IgamI(a, 0.75)
    45  	fhi := 0.25 - y
    46  
    47  	params := []float64{a, y}
    48  
    49  	// Also, after we generate a small interval by bisection above, false
    50  	// position will do a large step from an interval of width ~1e-4 to ~1e-14
    51  	// in one step (a=10, x=0.05, but similar for other values).
    52  	result, bestX, _, errEst := falsePosition(lo, hi, flo, fhi, 2*machEp, 2*machEp, 1e-2*a, gammaIncReg, params)
    53  	if result == fSolveMaxIterations && errEst > allowedATol+allowedRTol*math.Abs(bestX) {
    54  		bestX = math.NaN()
    55  	}
    56  
    57  	return bestX
    58  }