gonum.org/v1/gonum@v0.14.0/mathext/internal/cephes/igami.go (about)

     1  // Derived from SciPy's special/cephes/igami.c
     2  // https://github.com/scipy/scipy/blob/master/scipy/special/cephes/igami.c
     3  // Made freely available by Stephen L. Moshier without support or guarantee.
     4  
     5  // Use of this source code is governed by a BSD-style
     6  // license that can be found in the LICENSE file.
     7  // Copyright ©1984, ©1987, ©1995 by Stephen L. Moshier
     8  // Portions Copyright ©2017 The Gonum Authors. All rights reserved.
     9  
    10  package cephes
    11  
    12  import "math"
    13  
    14  // IgamI computes the inverse of the incomplete Gamma function. That is, it
    15  // returns the x such that:
    16  //
    17  //	IgamC(a, x) = p
    18  //
    19  // The input argument a must be positive and p must be between 0 and 1
    20  // inclusive or IgamI will panic. IgamI should return a positive number, but
    21  // can return 0 even with non-zero y due to underflow.
    22  func IgamI(a, p float64) float64 {
    23  	// Bound the solution
    24  	x0 := math.MaxFloat64
    25  	yl := 0.0
    26  	x1 := 0.0
    27  	yh := 1.0
    28  	dithresh := 5.0 * machEp
    29  
    30  	if p < 0 || p > 1 || a <= 0 {
    31  		panic(paramOutOfBounds)
    32  	}
    33  
    34  	if p == 0 {
    35  		return math.Inf(1)
    36  	}
    37  
    38  	if p == 1 {
    39  		return 0.0
    40  	}
    41  
    42  	// Starting with the approximate value
    43  	//  x = a y^3
    44  	// where
    45  	//  y = 1 - d - ndtri(p) sqrt(d)
    46  	// and
    47  	//  d = 1/9a
    48  	// the routine performs up to 10 Newton iterations to find the root of
    49  	//  IgamC(a, x) - p = 0
    50  	d := 1.0 / (9.0 * a)
    51  	y := 1.0 - d - Ndtri(p)*math.Sqrt(d)
    52  	x := a * y * y * y
    53  
    54  	lgm := lgam(a)
    55  
    56  	for i := 0; i < 10; i++ {
    57  		if x > x0 || x < x1 {
    58  			break
    59  		}
    60  
    61  		y = IgamC(a, x)
    62  
    63  		if y < yl || y > yh {
    64  			break
    65  		}
    66  
    67  		if y < p {
    68  			x0 = x
    69  			yl = y
    70  		} else {
    71  			x1 = x
    72  			yh = y
    73  		}
    74  
    75  		// Compute the derivative of the function at this point
    76  		d = (a-1)*math.Log(x) - x - lgm
    77  		if d < -maxLog {
    78  			break
    79  		}
    80  		d = -math.Exp(d)
    81  
    82  		// Compute the step to the next approximation of x
    83  		d = (y - p) / d
    84  		if math.Abs(d/x) < machEp {
    85  			return x
    86  		}
    87  		x = x - d
    88  	}
    89  
    90  	d = 0.0625
    91  	if x0 == math.MaxFloat64 {
    92  		if x <= 0 {
    93  			x = 1
    94  		}
    95  		for x0 == math.MaxFloat64 {
    96  			x = (1 + d) * x
    97  			y = IgamC(a, x)
    98  			if y < p {
    99  				x0 = x
   100  				yl = y
   101  				break
   102  			}
   103  			d = d + d
   104  		}
   105  	}
   106  
   107  	d = 0.5
   108  	dir := 0
   109  	for i := 0; i < 400; i++ {
   110  		x = x1 + d*(x0-x1)
   111  		y = IgamC(a, x)
   112  
   113  		lgm = (x0 - x1) / (x1 + x0)
   114  		if math.Abs(lgm) < dithresh {
   115  			break
   116  		}
   117  
   118  		lgm = (y - p) / p
   119  		if math.Abs(lgm) < dithresh {
   120  			break
   121  		}
   122  
   123  		if x <= 0 {
   124  			break
   125  		}
   126  
   127  		if y >= p {
   128  			x1 = x
   129  			yh = y
   130  			if dir < 0 {
   131  				dir = 0
   132  				d = 0.5
   133  			} else if dir > 1 {
   134  				d = 0.5*d + 0.5
   135  			} else {
   136  				d = (p - yl) / (yh - yl)
   137  			}
   138  			dir++
   139  		} else {
   140  			x0 = x
   141  			yl = y
   142  			if dir > 0 {
   143  				dir = 0
   144  				d = 0.5
   145  			} else if dir < -1 {
   146  				d = 0.5 * d
   147  			} else {
   148  				d = (p - yl) / (yh - yl)
   149  			}
   150  			dir--
   151  		}
   152  	}
   153  
   154  	return x
   155  }