github.com/gopherd/gonum@v0.0.4/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  //  IgamC(a, x) = p
    17  // The input argument a must be positive and p must be between 0 and 1
    18  // inclusive or IgamI will panic. IgamI should return a positive number, but
    19  // can return 0 even with non-zero y due to underflow.
    20  func IgamI(a, p float64) float64 {
    21  	// Bound the solution
    22  	x0 := math.MaxFloat64
    23  	yl := 0.0
    24  	x1 := 0.0
    25  	yh := 1.0
    26  	dithresh := 5.0 * machEp
    27  
    28  	if p < 0 || p > 1 || a <= 0 {
    29  		panic(paramOutOfBounds)
    30  	}
    31  
    32  	if p == 0 {
    33  		return math.Inf(1)
    34  	}
    35  
    36  	if p == 1 {
    37  		return 0.0
    38  	}
    39  
    40  	// Starting with the approximate value
    41  	//  x = a y^3
    42  	// where
    43  	//  y = 1 - d - ndtri(p) sqrt(d)
    44  	// and
    45  	//  d = 1/9a
    46  	// the routine performs up to 10 Newton iterations to find the root of
    47  	//  IgamC(a, x) - p = 0
    48  	d := 1.0 / (9.0 * a)
    49  	y := 1.0 - d - Ndtri(p)*math.Sqrt(d)
    50  	x := a * y * y * y
    51  
    52  	lgm := lgam(a)
    53  
    54  	for i := 0; i < 10; i++ {
    55  		if x > x0 || x < x1 {
    56  			break
    57  		}
    58  
    59  		y = IgamC(a, x)
    60  
    61  		if y < yl || y > yh {
    62  			break
    63  		}
    64  
    65  		if y < p {
    66  			x0 = x
    67  			yl = y
    68  		} else {
    69  			x1 = x
    70  			yh = y
    71  		}
    72  
    73  		// Compute the derivative of the function at this point
    74  		d = (a-1)*math.Log(x) - x - lgm
    75  		if d < -maxLog {
    76  			break
    77  		}
    78  		d = -math.Exp(d)
    79  
    80  		// Compute the step to the next approximation of x
    81  		d = (y - p) / d
    82  		if math.Abs(d/x) < machEp {
    83  			return x
    84  		}
    85  		x = x - d
    86  	}
    87  
    88  	d = 0.0625
    89  	if x0 == math.MaxFloat64 {
    90  		if x <= 0 {
    91  			x = 1
    92  		}
    93  		for x0 == math.MaxFloat64 {
    94  			x = (1 + d) * x
    95  			y = IgamC(a, x)
    96  			if y < p {
    97  				x0 = x
    98  				yl = y
    99  				break
   100  			}
   101  			d = d + d
   102  		}
   103  	}
   104  
   105  	d = 0.5
   106  	dir := 0
   107  	for i := 0; i < 400; i++ {
   108  		x = x1 + d*(x0-x1)
   109  		y = IgamC(a, x)
   110  
   111  		lgm = (x0 - x1) / (x1 + x0)
   112  		if math.Abs(lgm) < dithresh {
   113  			break
   114  		}
   115  
   116  		lgm = (y - p) / p
   117  		if math.Abs(lgm) < dithresh {
   118  			break
   119  		}
   120  
   121  		if x <= 0 {
   122  			break
   123  		}
   124  
   125  		if y >= p {
   126  			x1 = x
   127  			yh = y
   128  			if dir < 0 {
   129  				dir = 0
   130  				d = 0.5
   131  			} else if dir > 1 {
   132  				d = 0.5*d + 0.5
   133  			} else {
   134  				d = (p - yl) / (yh - yl)
   135  			}
   136  			dir++
   137  		} else {
   138  			x0 = x
   139  			yl = y
   140  			if dir > 0 {
   141  				dir = 0
   142  				d = 0.5
   143  			} else if dir < -1 {
   144  				d = 0.5 * d
   145  			} else {
   146  				d = (p - yl) / (yh - yl)
   147  			}
   148  			dir--
   149  		}
   150  	}
   151  
   152  	return x
   153  }