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  }