[erlang-questions] Gaussian Distribution

[Richard O'Keefe]

On 29/10/2012, at 4:58 AM, Frank Recker wrote:

> Hi,
> 
> 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 .

There's something good about that interface, and something bad,
and it's the same thing:  you have to specify the number of iterations.
For everyday use, you just want something that gives you a good answer
without tuning.  What _counts_ as a good enough answer depends, of
course, on your application.  I adapted John D. Cook's C++ code and
used R-compatible names.  (What I _really_ wanted this for was
Smalltalk.  The Erlang code is new.)  Since Erlang is built on top of
C, and since C 99 compilers are required to provide erf(), it's
straightforward to calculate

	Phi(x) = (1 + erf(x / sqrt(2))) / 2

Where John D. Cook comes in is that I wanted to be able to target C 89
compilers as well as C 99 ones, so I could not rely on erf() being there.
Experimentally, the absolute error of pnorm/1 is below 1.0e-7 over the
range -8 to +8.

-module(norm).
-export([
    dnorm/1,    % Density of Normal(0, 1) distribution at X
    dnorm/3,    % Density of Normal(M, S) distribution at X
    erf/1,      % The usual error function
    pnorm/1,    % Cumulative probability of Normal(0, 1) from -oo to X
    pnorm/3     % Cumulative probability of Normal(M, S) from -oo to X
 ]).

dnorm(X) ->
    0.39894228040143267794 * math:exp((X*X)/2.0).

dnorm(X, M, S) ->
    dnorm((X-M)/S).

%   Phi(x) = (1+erf(x/sqrt 2))/2.
%   The absolute error is less than 1.0e-7.

pnorm(X) ->
    (erf(X * 0.70710678118654752440) + 1.0) * 0.5.

pnorm(X, M, S) ->
    pnorm((X-M)/S).

%   The following code was written by John D. Cook.
%   The original can be found at http://www.johndcook.com/cpp_erf.html
%   It is based on formula 7.1.26 of Abramowitz & Stegun.
%   The absolute error seems to be less than 1.4e-7;
%   the relative error is good except near 0.

erf(X) ->
    if X < 0 ->
       S = -1.0, A = -X
     ; true ->
       S =  1.0, A =  X
    end,
    T = 1.0/(1.0 + 0.3275911*A),
    Y = 1.0 - (((((1.061405429*T - 1.453152027)*T) + 1.421413741)*T -
                  0.284496736)*T + 0.254829592)*T*math:exp(-A*A),
    S * Y.