Editing Normal distribution (section)

== Computational methods ==
=== Generating values from normal distribution ===
[[File:Planche de Galton.jpg|thumb|250px|right|The [[bean machine]], a device invented by [[Francis Galton]], can be called the first generator of normal random variables. This machine consists of a vertical board with interleaved rows of pins. Small balls are dropped from the top and then bounce randomly left or right as they hit the pins. The balls are collected into bins at the bottom and settle down into a pattern resembling the Gaussian curve.]]

In computer simulations, especially in applications of the [[Monte-Carlo method]], it is often desirable to generate values that are normally distributed. The algorithms listed below all generate the standard normal deviates, since a {{math|''N''(''μ'', ''σ''<sup>2</sup>)}} can be generated as {{math|1=''X'' = ''μ'' + ''σZ''}}, where ''Z'' is standard normal. All these algorithms rely on the availability of a [[random number generator]] ''U'' capable of producing [[Uniform distribution (continuous)|uniform]] random variates.
* The most straightforward method is based on the [[probability integral transform]] property: if ''U'' is distributed uniformly on (0,1), then Φ<sup>−1</sup>(''U'') will have the standard normal distribution. The drawback of this method is that it relies on calculation of the [[probit function]] Φ<sup>−1</sup>, which cannot be done analytically. Some approximate methods are described in {{harvtxt |Hart |1968 }} and in the [[error function|erf]] article. Wichura gives a fast algorithm for computing this function to 16 decimal places,<ref>{{cite journal|last=Wichura |first=Michael J.|year=1988|title=Algorithm AS241: The Percentage Points of the Normal Distribution|journal=Applied Statistics |volume=37|pages=477–84|doi=10.2307/2347330|jstor=2347330|issue=3}}</ref> which is used by [[R programming language|R]] to compute random variates of the normal distribution.
* [[Irwin–Hall distribution#Approximating a Normal distribution|An easy-to-program approximate approach]] that relies on the [[central limit theorem]] is as follows: generate 12 uniform ''U''(0,1) deviates, add them all up, and subtract 6 – the resulting random variable will have approximately standard normal distribution. In truth, the distribution will be [[Irwin–Hall distribution|Irwin–Hall]], which is a 12-section eleventh-order polynomial approximation to the normal distribution. This random deviate will have a limited range of (−6,&nbsp;6).<ref>{{harvtxt |Johnson |Kotz |Balakrishnan |1995 |loc=Equation (26.48) }}</ref> Note that in a true normal distribution, only 0.00034% of all samples will fall outside ±6σ.
* The [[Box–Muller method]] uses two independent random numbers ''U'' and ''V'' distributed [[uniform distribution (continuous)|uniformly]] on (0,1). Then the two random variables ''X'' and ''Y'' <math display=block>
    X = \sqrt{- 2 \ln U} \, \cos(2 \pi V) , \qquad
    Y = \sqrt{- 2 \ln U} \, \sin(2 \pi V) .
</math> will both have the standard normal distribution, and will be [[independence (probability theory)|independent]]. This formulation arises because for a [[bivariate normal]] random vector (''X'', ''Y'') the squared norm {{math|''X''<sup>2</sup> + ''Y''<sup>2</sup>}} will have the [[chi-squared distribution]] with two degrees of freedom, which is an easily generated [[exponential random variable]] corresponding to the quantity −2&nbsp;ln(''U'') in these equations; and the angle is distributed uniformly around the circle, chosen by the random variable ''V''.
* The [[Marsaglia polar method]] is a modification of the Box–Muller method which does not require computation of the sine and cosine functions. In this method, ''U'' and ''V'' are drawn from the uniform (−1,1) distribution, and then {{math|1=''S'' = ''U''<sup>2</sup> + ''V''<sup>2</sup>}} is computed. If ''S'' is greater or equal to 1, then the method starts over, otherwise the two quantities <math display=block>X = U\sqrt{\frac{-2\ln S}{S}}, \qquad Y = V\sqrt{\frac{-2\ln S}{S}}</math> are returned. Again, ''X'' and ''Y'' are independent, standard normal random variables.
* The Ratio method<ref>{{harvtxt |Kinderman |Monahan |1977 }}</ref> is a rejection method. The algorithm proceeds as follows:
** Generate two independent uniform deviates ''U'' and ''V'';
** Compute ''X'' = {{sqrt|8/''e''}} (''V'' − 0.5)/''U'';
** Optional: if ''X''<sup>2</sup> ≤ 5 − 4''e''<sup>1/4</sup>''U'' then accept ''X'' and terminate algorithm;
** Optional: if ''X''<sup>2</sup> ≥ 4''e''<sup>−1.35</sup>/''U'' + 1.4 then reject ''X'' and start over from step 1;
** If ''X''<sup>2</sup> ≤ −4 ln''U'' then accept ''X'', otherwise start over the algorithm.
*: The two optional steps allow the evaluation of the logarithm in the last step to be avoided in most cases. These steps can be greatly improved<ref>{{harvtxt|Leva|1992}}</ref> so that the logarithm is rarely evaluated.
* The [[ziggurat algorithm]]<ref>{{harvtxt |Marsaglia |Tsang |2000 }}</ref> is faster than the Box–Muller transform and still exact. In about 97% of all cases it uses only two random numbers, one random integer and one random uniform, one multiplication and an if-test. Only in 3% of the cases, where the combination of those two falls outside the "core of the ziggurat" (a kind of rejection sampling using logarithms), do exponentials and more uniform random numbers have to be employed.
* Integer arithmetic can be used to sample from the standard normal distribution.<ref>{{harvtxt|Karney|2016}}</ref><ref>{{harvtxt|Du|Fan|Wei|2022}}</ref> This method is exact in the sense that it satisfies the conditions of ''ideal approximation'';<ref>{{harvtxt|Monahan|1985|loc=section 2}}</ref> i.e., it is equivalent to sampling a real number from the standard normal distribution and rounding this to the nearest representable floating point number.
* There is also some investigation<ref>{{harvtxt |Wallace |1996}}</ref> into the connection between the fast [[Hadamard transform]] and the normal distribution, since the transform employs just addition and subtraction and by the central limit theorem random numbers from almost any distribution will be transformed into the normal distribution. In this regard a series of Hadamard transforms can be combined with random permutations to turn arbitrary data sets into a normally distributed data.

=== Numerical approximations for the normal cumulative distribution function and normal quantile function ===
The standard normal [[cumulative distribution function]] is widely used in scientific and statistical computing.

The values Φ(''x'') may be approximated very accurately by a variety of methods, such as [[numerical integration]], [[Taylor series]], [[asymptotic series]] and [[Gauss's continued fraction#Of Kummer's confluent hypergeometric function|continued fractions]]. Different approximations are used depending on the desired level of accuracy.
* {{harvtxt |Zelen |Severo |1964 }} give the approximation for Φ(''x'') for ''x'' > 0 with the absolute error {{math|{{abs|''ε''(''x'')}} < 7.5·10<sup>−8</sup>}} (algorithm [https://secure.math.ubc.ca/~cbm/aands/page_932.htm 26.2.17]): <math display=block>
    \Phi(x) = 1 - \varphi(x)\left(b_1 t + b_2 t^2 + b_3t^3 + b_4 t^4 + b_5 t^5\right) + \varepsilon(x), \qquad t = \frac{1}{1+b_0x},
</math> where ''ϕ''(''x'') is the standard normal probability density function, and ''b''<sub>0</sub> = 0.2316419, ''b''<sub>1</sub> = 0.319381530, ''b''<sub>2</sub> = −0.356563782, ''b''<sub>3</sub> = 1.781477937, ''b''<sub>4</sub> = −1.821255978, ''b''<sub>5</sub> = 1.330274429.
* {{harvtxt |Hart |1968 }} lists some dozens of approximations – by means of rational functions, with or without exponentials – for the {{mono|erfc()}} function. His algorithms vary in the degree of complexity and the resulting precision, with maximum absolute precision of 24 digits. An algorithm by {{harvtxt |West |2009 }} combines Hart's algorithm 5666 with a [[continued fraction]] approximation in the tail to provide a fast computation algorithm with a 16-digit precision.
* {{harvtxt |Cody |1969 }} after recalling Hart68 solution is not suited for erf, gives a solution for both erf and erfc, with maximal relative error bound, via [[rational function|Rational Chebyshev Approximation]].
* {{harvtxt |Marsaglia |2004 }} suggested a simple algorithm{{NoteTag|For example, this algorithm is given in the article [[Bc programming language#A translated C function|Bc programming language]].}} based on the Taylor series expansion <math display=block>
    \Phi(x) = \frac12 + \varphi(x)\left( x + \frac{x^3} 3 + \frac{x^5}{3 \cdot 5} + \frac{x^7}{3 \cdot 5 \cdot 7} + \frac{x^9}{3 \cdot 5 \cdot 7 \cdot 9} + \cdots \right)
</math> for calculating {{math|Φ(''x'')}} with arbitrary precision. The drawback of this algorithm is comparatively slow calculation time (for example it takes over 300 iterations to calculate the function with 16 digits of precision when {{math|1=''x'' = 10}}).
* The [[GNU Scientific Library]] calculates values of the standard normal cumulative distribution function using Hart's algorithms and approximations with [[Chebyshev polynomial]]s.
* {{harvtxt |Dia|2023 }} proposes the following approximation of <math display=inline>1-\Phi</math> with a maximum relative error less than <math display=inline>2^{-53}</math> <math display=inline>
\left(\approx 1.1 \times 10^{-16}\right)
</math> in absolute value: for <math display=inline>x \ge 0</math><math display=inline>
\begin{aligned}
    1-\Phi\left(x\right) & = \left(\frac{0.39894228040143268}{x+2.92678600515804815}\right)
    \left(\frac{x^2+8.42742300458043240 x+18.38871225773938487}{x^2+5.81582518933527391 x+8.97280659046817350} \right) \\
    & \left(\frac{x^2+7.30756258553673541 x+18.25323235347346525}{x^2+5.70347935898051437 x+10.27157061171363079}\right)
    \left(\frac{x^2+5.66479518878470765 x+18.61193318971775795}{x^2+5.51862483025707963 x+12.72323261907760928}\right) \\
    & \left( \frac{x^2+4.91396098895240075 x+24.14804072812762821}{x^2+5.26184239579604207 x+16.88639562007936908}\right)
    \left( \frac{x^2+3.83362947800146179 x+11.61511226260603247}{x^2+4.92081346632882033 x+24.12333774572479110}\right) e^{-\frac{x^2}{2}}
\end{aligned}
</math> and for <math display=inline>
    x<0
</math>,
<math display=block>
    1-\Phi\left(x\right) = 1-\left(1-\Phi\left(-x\right)\right)
</math>

Shore (1982) introduced simple approximations that may be incorporated in stochastic optimization models of engineering and operations research, like reliability engineering and inventory analysis. Denoting {{math|1=''p'' = Φ(''z'')}}, the simplest approximation for the quantile function is:
<math display=block>z = \Phi^{-1}(p)=5.5556\left[1- \left( \frac{1-p} p \right)^{0.1186}\right],\qquad p\ge 1/2</math>

This approximation delivers for ''z'' a maximum absolute error of 0.026 (for {{math|0.5 ≤ ''p'' ≤ 0.9999}}, corresponding to {{math|0 ≤ ''z'' ≤ 3.719}}). For {{math|''p'' < 1/2}} replace ''p'' by {{math|1 − ''p''}} and change sign. Another approximation, somewhat less accurate, is the single-parameter approximation:
<math display=block> z=-0.4115\left\{ \frac{1-p} p + \log \left[ \frac{1-p} p \right] - 1 \right\}, \qquad p\ge 1/2</math>

The latter had served to derive a simple approximation for the loss integral of the normal distribution, defined by
<math display=block>\begin{align}
L(z) & =\int_z^\infty (u-z)\varphi(u) \, du=\int_z^\infty [1-\Phi (u)] \, du \\[5pt]
L(z) & \approx \begin{cases}
   0.4115\left(\dfrac p {1-p} \right) - z, & p<1/2, \\ \\
   0.4115\left( \dfrac {1-p} p \right), & p\ge 1/2.
\end{cases} \\[5pt]
\text{or, equivalently,} \\
L(z) & \approx \begin{cases}
   0.4115\left\{ 1-\log \left[ \frac p {1-p} \right] \right\}, & p < 1/2, \\ \\
   0.4115 \dfrac{1-p} p, & p\ge 1/2.
\end{cases}
\end{align}</math>

This approximation is particularly accurate for the right far-tail (maximum error of 10<sup>−3</sup> for z≥1.4). Highly accurate approximations for the cumulative distribution function, based on [[Response Modeling Methodology]] (RMM, Shore, 2011, 2012), are shown in Shore (2005).

Some more approximations can be found at: [[Error function#Approximation with elementary functions]]. In particular, small ''relative'' error on the whole domain for the cumulative distribution function {{tmath|\Phi}} and the quantile function <math display=inline>\Phi^{-1}</math> as well, is achieved via an explicitly invertible formula by Sergei Winitzki in 2008.