Xref: utzoo sci.math.stat:2128 comp.lang.pascal:5987
Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!usc!alhena.usc.edu!ajayshah
From: ajayshah@alhena.usc.edu (Ajay Shah)
Newsgroups: sci.math.stat,comp.lang.pascal
Subject: TP code (Was Re: Any normally-distributed random number functions?)
Keywords: random generator normal density algorithm
Message-ID: <31641@usc>
Date: 4 Apr 91 22:19:34 GMT
References: <1991Apr04.183739.4876@group1.UUCP>
Sender: news@usc
Followup-To: sci.math.stat
Organization: University of Southern California, Los Angeles, CA
Lines: 275
Nntp-Posting-Host: alhena.usc.edu


{$N+} {$E+}

UNIT distrns;
     {All kind of code connected with various distributions.}

INTERFACE

{NORMAL DISTRIBUTION}
{3 kernel functions:}
function StdNorm(x:double):double; {standard normal density function}
function CNorm(x:double):double; {cumulative standard normal df}
function BivNormDens(X, Y, rho:double):double; {bivariate normal density}

{oft-used gloss}
function lndnorm(x:double):double; {compute ln(stdnormdens at x), in 3 flops!}
function NormalDensity(mu, sigma, x:double):double;
function NormArea(x1, x2:double):double; {area between x1 and x2 on N(0,1)}
function mi(x:double):double; {ratio of density upon area}

{random numbers:}
function StdRandNorm:double;
function RandNorm(myoo, sigma:double):double;
procedure StdBivNorRand(rho:double; var X1, X2:double);
          {returns draw from mean zero, variance 1, corr. coeff rho
           bivariate normal distribution.}

{FEW OTHERS}

function CauchyRand(var h, m:double):double;
function ExpRand(var a:double):double;
function Uniform(var a,b:double):double;

IMPLEMENTATION

const
     NORMMAX = 20.0; {place where we force normal density to zero,
                      similar restrictions on normal probability.}
     SQRT2 = 1.414213562373095049;
     SQRTPI = 1.772453850905516027;
     OneUponSqrt2Pi = 0.398942280401433;
     TWOPI = 6.2831853071795864769;
     LNSQRT2PI = -0.9189385332046727417803296;


{----------------------------------------------------------------------------}
function StdNorm(x:double):double; {standard normal density function}
begin
     if abs(x) > NORMMAX then StdNorm := 0.0
                         else StdNorm := exp(-0.5*x*x)*OneUponSqrt2Pi
end;

{----------------------------------------------------------------------------}
function NormalDensity(mu, sigma, x:double):double;
var tmp:double;
begin
     tmp := (x-mu)/sigma;
     NormalDensity := StdNorm(tmp)/sigma
end;


{----------------------------------------------------------------------------}
function CNorm(x:double):double;

{Objective :
           given x, returns integral from -infinity to x of the std.
           normal distribution.

 Precision:
           all the digits printed in Abramovitch and Stegun (14? 17?)

 Method :
           it knows 3 approximations, in the regions
              (I)   abs(x/root2) < 0.46
              (II)  0.46 < abs(x/root2) < 4
              (III) 4 < abs(x/root2)
           in each of these regions, it approximates CNorm by a polynomial-ish.

 Speed :
           On my 386/387 at 20 MHz, no cache, Nortons SI = 21.1,
           region (I):   0.2 milliseconds
           region (II):  0.47 milliseconds
           region (III): 0.43 milliseconds

 History :
           * Originally a CACM algorithm, forget which,
           * Translated from fortran into C by Luke Tierney (XLispStat),
           * Translated from C into Turbo Pascal by Ajay Shah, 14 Dec 1990.
 }

const
	P10 = 242.66795523053175; P11 = 21.979261618294152;
	P12 = 6.9963834886191355; P13 = -0.035609843701815385;
	Q10 = 215.05887586986120; Q11 = 91.164905404514901;
	Q12 = 15.082797630407787; Q13 = 1.0;

	P20 = 300.4592610201616005;  P21 = 451.9189537118729422;
	P22 = 339.3208167343436870;  P23 = 152.9892850469404039;
	P24 = 43.16222722205673530;  P25 = 7.211758250883093659;
	P26 = 0.5641955174789739711; P27 = -0.0000001368648573827167067;
	Q20 = 300.4592609569832933;  Q21 = 790.9509253278980272;
	Q22 = 931.3540948506096211;  Q23 = 638.9802644656311665;
	Q24 = 277.5854447439876434;  Q25 = 77.00015293522947295;
	Q26 = 12.78272731962942351;  Q27 = 1.0;

	P30 = -0.00299610707703542174; P31 = -0.0494730910623250734;
	P32 = -0.226956593539686930;   P33 = -0.278661308609647788;
	P34 = -0.0223192459734184686;  Q30 = 0.0106209230528467918;
	Q31 = 0.191308926107829841;    Q32 = 1.05167510706793207;
	Q33 = 1.98733201817135256;     Q34 = 1.0;
{These constants are totally specific to this function; hence they're
 not at the top of this file with sensible names.}

var
   sn:integer;
   R1, R2, R3, y, y2, y3, y4, y5, y6, y7, erf, erfc, z, z2, z3, z4 : double;

begin
     if (x < -NORMMAX) then
     begin
          CNorm := 0.0; exit
     end;
     if (x > NORMMAX) then
     begin
          CNorm := 1.0; exit
     end;

     y := x/SQRT2;
     if (y < 0.0) then begin
                            y := -y;
                            sn := -1;
                       end
                  else sn := 1;

     y2 := y * y;
     y4 := y2 * y2;
     y6 := y4 * y2;

     if (y < 0.46875) then
     begin
          R1 := P10 + (P11*y2) + (P12*y4) + (P13*y6);
          R2 := Q10 + (Q11*y2) + (Q12*y4) + (Q13*y6);
          erf := y * R1 / R2;
          if (sn = 1) then CNorm := 0.5 + (erf/2.0)
                      else CNorm := 0.5 - (erf/2.0);
          exit
     end;

     if (y < 4.0) then
     begin
          y3 := y2 * y;
          y5 := y4 * y;
          y7 := y6 * y;
          R1 := P20 + (P21*y) + (P22*y2) + (P23*y3) + (P24*y4) +
                (P25*y5) + (P26*y6) + (P27*y7);
          R2 := Q20 + (Q21*y) + (Q22*y2) + (Q23*y3) + (Q24*y4) +
                (Q25*y5) + (Q26*y6) + (Q27*y7);
          erfc := exp(-y2) * R1 / R2;
          if (sn = 1) then CNorm := 1.0 - (erfc/2.0)
                      else CNorm := erfc/2.0;
          exit
     end;

     z := y4;
     z2 := z * z;
     z3 := z2 * z;
     z4 := z2 * z2;
     R1 := P30 + (P31*z) + (P32*z2) + (P33*z3) + (P34*z4);
     R2 := Q30 + (Q31*z) + (Q32*z2) + (Q33*z3) + (Q34*z4);
     erfc := (exp(-y2)/y) * ( (1.0/SQRTPI) + (R1/(R2 * y2)) );
     if (sn = 1) then CNorm := 1.0 - (erfc/2.0)
                 else CNorm := erfc/2.0;
end;


{----------------------------------------------------------------------------}
function lndnorm(x:double):double; {compute ln(stdnormdens at x), in 3 flops!}
begin
     lndnorm :=  LNSQRT2PI - (sqr(x)/2.0)
end;


{----------------------------------------------------------------------------}
function NormArea(x1, x2:double):double; {area between x1 and x2 on N(0,1)}
begin {You just get junk if x1 > x2}
     NormArea := CNorm(x2) - CNorm(x1)
end;

{----------------------------------------------------------------------------}
function mi(x:double):double; {ratio of density upon area}
begin
     mi := StdNorm(x)/cnorm(x)
end;


{----------------------------------------------------------------------------}
function bivnormdens(X, Y, rho:double):double;
{(X,Y) is a point in x-y plane, rho = correlation coefficient.
 The function computes bivariate std. normal density at (x,y).}
var
   rho2, tmp, denom, power : double;
begin
      rho2 := sqr(rho);
      tmp := 1.0 - rho2;
      denom := TWOPI*sqrt(tmp);
      power := (X*(X - (2.0*rho*Y)) + sqr(Y))/tmp;
      bivnormdens := exp(-power/2.0) / denom
end;


{----------------------------------------------------------------------------}
function stdrandnorm:double;
{reference: page 953 of Abramowitch and Stegun}
begin
     stdrandnorm := sqrt(-2.0*ln(random))*cos(TWOPI*random)
end;


{----------------------------------------------------------------------------}
function randnorm(myoo, sigma:double):double;

begin
     randnorm := (sigma*stdrandnorm) + myoo
end;

{----------------------------------------------------------------------------}
function tan(x:double):double;
begin
     tan := sin(x)/cos(x)
end;

{----------------------------------------------------------------------------}
function CauchyRand(var h, m:double):double;
begin
     CauchyRand := h*(tan(pi*(random - 0.5))) + m
end;

{----------------------------------------------------------------------------}
function ExpRand(var a:double):double;
begin
     ExpRand := -ln(random)/a;
end;

{----------------------------------------------------------------------------}
function Uniform(var a,b:double):double;
begin
     Uniform := a + (random*(b-a))
end;

{----------------------------------------------------------------------------}
procedure StdBivNorRand(rho:double; var X1, X2:double);
{reference: page 953 of Abramowitch and Stegun}
var U1, U2, tmp1, tmp2 : double;
begin
     U1 := random; U2 := random;
        {generating uniformly distributed random numbers ain't cheap!}
     tmp1 := sqrt(-2.0*ln(U1));
     tmp2 := TWOPI*U2;
     X1 := tmp1*cos(tmp2);
     X2 := tmp1*sin(tmp2);
     {X1 and X2 are now independent, std. normal random numbers.
      That was cheaper than calling RandNorm twice.}
     X2 := (rho*X1) + (X2*sqrt(1.0-(rho*rho)))
     {X1 and X2 are now drawn from the required bivariate normal}
end;

begin {unit init}
     randomize
end.

-- 
_______________________________________________________________________________
Ajay Shah, (213)734-3930, ajayshah@usc.edu
                             The more things change, the more they stay insane.
_______________________________________________________________________________