[erlang-questions] Gaussian Distribution

[Jachym Holecek]

# 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:      pjacklam@REDACTED
%% 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).