# [erlang-questions] Gaussian Distribution

Jachym Holecek <>
Mon Oct 29 10:31:33 CET 2012

```# Frank Recker 2012-10-28:
> at work, I often need the values the cumulative distribution function of
> the Gaussian distribution. The code for this function in haskell, erlang
> and perl and the corresponding mathematical paper can be found at
> git://github.com/frecker/gaussian-distribution.git .

One more to the mix -- inverse CDF for normal distribution. Only mildly
tested.

BR,
-- Jachym

%% Function normal_quant/1 is and Erlang version of:
%%
%%   http://home.online.no/~pjacklam/notes/invnorm/impl/acklam/bc/ltqnorm.bc
%%
%% 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.
%%
%% The shell syntax is (or should be, anyway) POSIX bc.
%%
%% Author:      Peter John Acklam
%% Time-stamp:  2005-03-10 14:13:52 +01:00
%% E-mail:
%% WWW URL:     http://home.online.no/~pjacklam

%% Inverse CDF of standard normal distribution.
normal_quant(P) when is_float(P) ->
Q = P - 0.5,
A = quant_fabs(Q),
if A =< 0.47575 ->
%% Rational approximation for central region.
R = Q*Q,
((((( -39.69683028665376*R + 220.9460984245205)*R - 275.9285104469687)*R +
138.3577518672690)*R - 30.66479806614716)*R + 2.506628277459239)*Q /
((((( -54.47609879822406*R + 161.5858368580409)*R - 155.6989798598866)*R +
66.80131188771972)*R - 13.28068155288572)*R + 1);
A > 0.47575 ->
%% Rational approximation for tails.
if Q > 0 ->
-quant_tails(1 - P);
true ->
quant_tails(P)
end
end.

quant_tails(P) when is_float(P) ->
R = math:sqrt(-2 * math:log(P)),
(((((-0.007784894002430293*R - 0.3223964580411365)*R - 2.400758277161838)*R -
2.549732539343734)*R + 4.374664141464968) * R + 2.938163982698783) /
((((  0.007784695709041462*R + 0.3224671290700398)*R + 2.445134137142996)*R +
3.754408661907416)*R + 1).

quant_fabs(X) when is_float(X), X < 0.0 ->
-X;
quant_fabs(X) when is_float(X) ->
X.

%% Inverse CDF of generic normal distribution.
normal_quant(P, Mean, Sigma) when is_float(Mean), is_float(Sigma) ->
Mean + math:sqrt(Sigma)*normal_quant(P).

```