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

     1  // Copyright ©2016 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  /*
     6   * Cephes Math Library Release 2.1:  January, 1989
     7   * Copyright 1984, 1987, 1989 by Stephen L. Moshier
     8   * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
     9   */
    10  
    11  package cephes
    12  
    13  import "math"
    14  
    15  // TODO(btracey): There is currently an implementation of this functionality
    16  // in gonum/stat/distuv. Find out which implementation is better, and rectify
    17  // by having distuv call this, or moving this implementation into
    18  // gonum/mathext/internal/gonum.
    19  
    20  // math.Sqrt(2*pi)
    21  const s2pi = 2.50662827463100050242e0
    22  
    23  // approximation for 0 <= |y - 0.5| <= 3/8
    24  var P0 = [5]float64{
    25  	-5.99633501014107895267e1,
    26  	9.80010754185999661536e1,
    27  	-5.66762857469070293439e1,
    28  	1.39312609387279679503e1,
    29  	-1.23916583867381258016e0,
    30  }
    31  
    32  var Q0 = [8]float64{
    33  	/* 1.00000000000000000000E0, */
    34  	1.95448858338141759834e0,
    35  	4.67627912898881538453e0,
    36  	8.63602421390890590575e1,
    37  	-2.25462687854119370527e2,
    38  	2.00260212380060660359e2,
    39  	-8.20372256168333339912e1,
    40  	1.59056225126211695515e1,
    41  	-1.18331621121330003142e0,
    42  }
    43  
    44  // Approximation for interval z = math.Sqrt(-2 log y ) between 2 and 8
    45  // i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
    46  var P1 = [9]float64{
    47  	4.05544892305962419923e0,
    48  	3.15251094599893866154e1,
    49  	5.71628192246421288162e1,
    50  	4.40805073893200834700e1,
    51  	1.46849561928858024014e1,
    52  	2.18663306850790267539e0,
    53  	-1.40256079171354495875e-1,
    54  	-3.50424626827848203418e-2,
    55  	-8.57456785154685413611e-4,
    56  }
    57  
    58  var Q1 = [8]float64{
    59  	/*  1.00000000000000000000E0, */
    60  	1.57799883256466749731e1,
    61  	4.53907635128879210584e1,
    62  	4.13172038254672030440e1,
    63  	1.50425385692907503408e1,
    64  	2.50464946208309415979e0,
    65  	-1.42182922854787788574e-1,
    66  	-3.80806407691578277194e-2,
    67  	-9.33259480895457427372e-4,
    68  }
    69  
    70  // Approximation for interval z = math.Sqrt(-2 log y ) between 8 and 64
    71  // i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
    72  var P2 = [9]float64{
    73  	3.23774891776946035970e0,
    74  	6.91522889068984211695e0,
    75  	3.93881025292474443415e0,
    76  	1.33303460815807542389e0,
    77  	2.01485389549179081538e-1,
    78  	1.23716634817820021358e-2,
    79  	3.01581553508235416007e-4,
    80  	2.65806974686737550832e-6,
    81  	6.23974539184983293730e-9,
    82  }
    83  
    84  var Q2 = [8]float64{
    85  	/*  1.00000000000000000000E0, */
    86  	6.02427039364742014255e0,
    87  	3.67983563856160859403e0,
    88  	1.37702099489081330271e0,
    89  	2.16236993594496635890e-1,
    90  	1.34204006088543189037e-2,
    91  	3.28014464682127739104e-4,
    92  	2.89247864745380683936e-6,
    93  	6.79019408009981274425e-9,
    94  }
    95  
    96  // Ndtri returns the argument, x, for which the area under the
    97  // Gaussian probability density function (integrated from
    98  // minus infinity to x) is equal to y.
    99  func Ndtri(y0 float64) float64 {
   100  	// For small arguments 0 < y < exp(-2), the program computes
   101  	// z = math.Sqrt( -2.0 * math.Log(y) );  then the approximation is
   102  	// x = z - math.Log(z)/z  - (1/z) P(1/z) / Q(1/z).
   103  	// There are two rational functions P/Q, one for 0 < y < exp(-32)
   104  	// and the other for y up to exp(-2).  For larger arguments,
   105  	// w = y - 0.5, and  x/math.Sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
   106  	var x, y, z, y2, x0, x1 float64
   107  	var code int
   108  
   109  	if y0 <= 0.0 {
   110  		if y0 < 0 {
   111  			panic(paramOutOfBounds)
   112  		}
   113  		return math.Inf(-1)
   114  	}
   115  	if y0 >= 1.0 {
   116  		if y0 > 1 {
   117  			panic(paramOutOfBounds)
   118  		}
   119  		return math.Inf(1)
   120  	}
   121  	code = 1
   122  	y = y0
   123  	if y > (1.0 - 0.13533528323661269189) { /* 0.135... = exp(-2) */
   124  		y = 1.0 - y
   125  		code = 0
   126  	}
   127  
   128  	if y > 0.13533528323661269189 {
   129  		y = y - 0.5
   130  		y2 = y * y
   131  		x = y + y*(y2*polevl(y2, P0[:], 4)/p1evl(y2, Q0[:], 8))
   132  		x = x * s2pi
   133  		return (x)
   134  	}
   135  
   136  	x = math.Sqrt(-2.0 * math.Log(y))
   137  	x0 = x - math.Log(x)/x
   138  
   139  	z = 1.0 / x
   140  	if x < 8.0 { /* y > exp(-32) = 1.2664165549e-14 */
   141  		x1 = z * polevl(z, P1[:], 8) / p1evl(z, Q1[:], 8)
   142  	} else {
   143  		x1 = z * polevl(z, P2[:], 8) / p1evl(z, Q2[:], 8)
   144  	}
   145  	x = x0 - x1
   146  	if code != 0 {
   147  		x = -x
   148  	}
   149  	return (x)
   150  }