github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/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  	"github.com/jingcheng-WU/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  //  Igam(a, x) = y
    28  // The input argument a must be positive and y must be between 0 and 1
    29  // inclusive or gammaIncRegInv will panic. gammaIncRegInv should return a
    30  // positive number, but can return NaN if there is a failure to converge.
    31  func gammaIncRegInv(a, y float64) float64 {
    32  	// For y not small, we just use
    33  	//  IgamI(a, 1-y)
    34  	// (inverse of the complemented incomplete Gamma integral). For y small,
    35  	// however, 1-y is about 1, and we lose digits.
    36  	if a <= 0 || y <= 0 || y >= 0.25 {
    37  		return cephes.IgamI(a, 1-y)
    38  	}
    39  
    40  	lo := 0.0
    41  	flo := -y
    42  	hi := cephes.IgamI(a, 0.75)
    43  	fhi := 0.25 - y
    44  
    45  	params := []float64{a, y}
    46  
    47  	// Also, after we generate a small interval by bisection above, false
    48  	// position will do a large step from an interval of width ~1e-4 to ~1e-14
    49  	// in one step (a=10, x=0.05, but similar for other values).
    50  	result, bestX, _, errEst := falsePosition(lo, hi, flo, fhi, 2*machEp, 2*machEp, 1e-2*a, gammaIncReg, params)
    51  	if result == fSolveMaxIterations && errEst > allowedATol+allowedRTol*math.Abs(bestX) {
    52  		bestX = math.NaN()
    53  	}
    54  
    55  	return bestX
    56  }