github.com/jgbaldwinbrown/perf@v0.1.1/pkg/stats/normaldist.go (about) 1 // Copyright 2015 The Go 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 package stats 6 7 import ( 8 "math" 9 "math/rand" 10 ) 11 12 // NormalDist is a normal (Gaussian) distribution with mean Mu and 13 // standard deviation Sigma. 14 type NormalDist struct { 15 Mu, Sigma float64 16 } 17 18 // StdNormal is the standard normal distribution (Mu = 0, Sigma = 1) 19 var StdNormal = NormalDist{0, 1} 20 21 // 1/sqrt(2 * pi) 22 const invSqrt2Pi = 0.39894228040143267793994605993438186847585863116493465766592583 23 24 func (n NormalDist) PDF(x float64) float64 { 25 z := x - n.Mu 26 return math.Exp(-z*z/(2*n.Sigma*n.Sigma)) * invSqrt2Pi / n.Sigma 27 } 28 29 func (n NormalDist) pdfEach(xs []float64) []float64 { 30 res := make([]float64, len(xs)) 31 if n.Mu == 0 && n.Sigma == 1 { 32 // Standard normal fast path 33 for i, x := range xs { 34 res[i] = math.Exp(-x*x/2) * invSqrt2Pi 35 } 36 } else { 37 a := -1 / (2 * n.Sigma * n.Sigma) 38 b := invSqrt2Pi / n.Sigma 39 for i, x := range xs { 40 z := x - n.Mu 41 res[i] = math.Exp(z*z*a) * b 42 } 43 } 44 return res 45 } 46 47 func (n NormalDist) CDF(x float64) float64 { 48 return math.Erfc(-(x-n.Mu)/(n.Sigma*math.Sqrt2)) / 2 49 } 50 51 func (n NormalDist) cdfEach(xs []float64) []float64 { 52 res := make([]float64, len(xs)) 53 a := 1 / (n.Sigma * math.Sqrt2) 54 for i, x := range xs { 55 res[i] = math.Erfc(-(x-n.Mu)*a) / 2 56 } 57 return res 58 } 59 60 func (n NormalDist) InvCDF(p float64) (x float64) { 61 // This is based on Peter John Acklam's inverse normal CDF 62 // algorithm: http://home.online.no/~pjacklam/notes/invnorm/ 63 const ( 64 a1 = -3.969683028665376e+01 65 a2 = 2.209460984245205e+02 66 a3 = -2.759285104469687e+02 67 a4 = 1.383577518672690e+02 68 a5 = -3.066479806614716e+01 69 a6 = 2.506628277459239e+00 70 71 b1 = -5.447609879822406e+01 72 b2 = 1.615858368580409e+02 73 b3 = -1.556989798598866e+02 74 b4 = 6.680131188771972e+01 75 b5 = -1.328068155288572e+01 76 77 c1 = -7.784894002430293e-03 78 c2 = -3.223964580411365e-01 79 c3 = -2.400758277161838e+00 80 c4 = -2.549732539343734e+00 81 c5 = 4.374664141464968e+00 82 c6 = 2.938163982698783e+00 83 84 d1 = 7.784695709041462e-03 85 d2 = 3.224671290700398e-01 86 d3 = 2.445134137142996e+00 87 d4 = 3.754408661907416e+00 88 89 plow = 0.02425 90 phigh = 1 - plow 91 ) 92 93 if p < 0 || p > 1 { 94 return nan 95 } else if p == 0 { 96 return -inf 97 } else if p == 1 { 98 return inf 99 } 100 101 if p < plow { 102 // Rational approximation for lower region. 103 q := math.Sqrt(-2 * math.Log(p)) 104 x = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q + c6) / 105 ((((d1*q+d2)*q+d3)*q+d4)*q + 1) 106 } else if phigh < p { 107 // Rational approximation for upper region. 108 q := math.Sqrt(-2 * math.Log(1-p)) 109 x = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q + c6) / 110 ((((d1*q+d2)*q+d3)*q+d4)*q + 1) 111 } else { 112 // Rational approximation for central region. 113 q := p - 0.5 114 r := q * q 115 x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r + a6) * q / 116 (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r + 1) 117 } 118 119 // Refine approximation. 120 e := 0.5*math.Erfc(-x/math.Sqrt2) - p 121 u := e * math.Sqrt(2*math.Pi) * math.Exp(x*x/2) 122 x = x - u/(1+x*u/2) 123 124 // Adjust from standard normal. 125 return x*n.Sigma + n.Mu 126 } 127 128 func (n NormalDist) Rand(r *rand.Rand) float64 { 129 var x float64 130 if r == nil { 131 x = rand.NormFloat64() 132 } else { 133 x = r.NormFloat64() 134 } 135 return x*n.Sigma + n.Mu 136 } 137 138 func (n NormalDist) Bounds() (float64, float64) { 139 const stddevs = 3 140 return n.Mu - stddevs*n.Sigma, n.Mu + stddevs*n.Sigma 141 }