-
Notifications
You must be signed in to change notification settings - Fork 146
/
ltqnorm.go
88 lines (82 loc) · 2.42 KB
/
ltqnorm.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
// Copyright 2016 The Gosl Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package rnd
import (
"math"
"github.com/cpmech/gosl/chk"
)
/* Lower tail quantile for standard normal distribution function.
*
* This function returns an approximation of the inverse cumulative
* standard normal distribution function. I.e., given P, it returns
* an approximation to the X satisfying P = Pr{Z <= X} where Z is a
* random variable from the standard normal distribution.
*
* The algorithm uses a minimax approximation by rational functions
* and the result has a relative error whose absolute value is less
* than 1.15e-9.
*
* Author: Peter John Acklam
* Time-stamp: 2002-06-09 18:45:44 +0200
* E-mail: [email protected]
* WWW URL: https://www.math.uio.no/~jacklam
*
* C implementation adapted from Peter's Perl version
*/
// coefficients in rational approximations
var (
a = []float64{
-3.969683028665376e+01,
2.209460984245205e+02,
-2.759285104469687e+02,
1.383577518672690e+02,
-3.066479806614716e+01,
2.506628277459239e+00,
}
b = []float64{
-5.447609879822406e+01,
1.615858368580409e+02,
-1.556989798598866e+02,
6.680131188771972e+01,
-1.328068155288572e+01,
}
c = []float64{
-7.784894002430293e-03,
-3.223964580411365e-01,
-2.400758277161838e+00,
-2.549732539343734e+00,
4.374664141464968e+00,
2.938163982698783e+00,
}
d = []float64{
7.784695709041462e-03,
3.224671290700398e-01,
2.445134137142996e+00,
3.754408661907416e+00,
}
)
func ltqnorm(p float64) float64 {
low := 0.02425
hight := 0.97575
var q, r float64
if p < 0 || p > 1 {
chk.Panic("input value (%g) is outside allowed range [0, 1]", p)
} else if p == 0 {
chk.Panic("p==0 => -infinity")
} else if p == 1 {
chk.Panic("p==1 => infinity")
} else if p < low {
// rational approximation for lower region
q = math.Sqrt(-2 * math.Log(p))
return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q + c[5]) / ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q + 1)
} else if p > hight {
// rational approximation for upper region
q = math.Sqrt(-2 * math.Log(1-p))
return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q + c[5]) / ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q + 1)
}
// rational approximation for central region
q = p - 0.5
r = q * q
return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r + a[5]) * q / (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r + 1)
}