ajayshah@alhena.usc.edu (Ajay Shah) (04/05/91)
{$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.
_______________________________________________________________________________