Category : Science and Education
Archive   : PROB.ZIP
Filename : PROB.PAS

PROGRAM Prob;
{ Given the mean, standard deviation, and upper and lower limits of a Gaussian }
{ distribution, find the associated probability. If executed from DOS as }
{ "prob", 6 digits accuracy are given. Other accuracies are found by passing }
{ a tolerance such as "prob 1.0e-15". Written by: Tom Hardy 9 March 1991 }
{$N+}
{$E+}
VAR
Quit:BOOLEAN;
Code:INTEGER;
Mean,sd,a,b,Error:EXTENDED;
ErrorString:STRING;

FUNCTION Probability(Lower,Upper,Error:EXTENDED):EXTENDED;
VAR
Temp:EXTENDED;

FUNCTION ERF(a:EXTENDED):EXTENDED;
{ Compute the Error Function integral from 0 to a, where a>=0 }
{ by integrating the Error Function }
VAR
i,n:LONGINT;
Sum,LastSum,h,Delta,x:EXTENDED;
BEGIN { ERF }
n:=100;
LastSum:=-1.0;
REPEAT
n:=2*n;
h:=a/(2*n);
Sum:=1.0;
FOR i:=1 TO 2*n-1 DO
BEGIN
x:=i*h;
IF ODD(i) THEN
Sum:=Sum+4*EXP(-SQR(x)/2)
ELSE
Sum:=Sum+2*EXP(-SQR(x)/2)
END;
Sum:=(h/3)*(Sum+EXP(-SQR(a)/2));
Delta:=ABS(Sum-LastSum);
LastSum:=Sum
UNTIL(Delta<=Error);
ERF:=Sum/SQRT(2*Pi)
END; { ERF }

BEGIN { Probability }
IF (Lower>=Upper) THEN
Probability:=0.0
ELSE
BEGIN
IF (Lower<-10.0) THEN
Lower:=-10.0;
IF (Upper>10.0) THEN
Upper:=10.0;
IF (Lower=0.0) THEN
Temp:=ERF(Upper)
ELSE
BEGIN
IF (Upper=0.0) THEN
Temp:=ERF(-Lower)
ELSE
BEGIN
IF (Lower>0.0) THEN
Temp:=ERF(Upper)-ERF(Lower)
ELSE
BEGIN
IF (Upper<=0.0) THEN
Temp:=ERF(-Lower)-ERF(-Upper)
ELSE
Temp:=ERF(-Lower)+ERF(Upper)
END
END
END;
IF (Temp>Error) THEN
Probability:=Temp
ELSE
Probability:=0.0
END
END; { Probability }

PROCEDURE GetReal(Str:STRING;VAR x:EXTENDED;VAR Quit:BOOLEAN);
VAR
Valid:BOOLEAN;
Code:INTEGER;
s:STRING;
BEGIN { GetReal }
REPEAT
WRITE(Str);
READLN(s);
Quit:=(s='');
Val(s,x,Code);
Valid:=(Code=0)
UNTIL (Valid OR Quit)
END; { GetReal }

PROCEDURE Getsd(VAR sd:EXTENDED;VAR Quit:BOOLEAN);
VAR
Error:BOOLEAN;
BEGIN { Getsd }
REPEAT
GetReal('Enter distribution standard deviation : ',sd,Quit);
Error:=(sd<=1.0e-3);
IF (sd<0.0) THEN
WRITELN('Standard deviation must be positive')
ELSE
IF (sd<=1.0E-3) THEN
WRITELN('Standard deviation too small')
UNTIL (Quit OR (NOT Error))
END; { Getsd }

PROCEDURE GetDataAndResults(Error:EXTENDED);
VAR
Digits:INTEGER;
BEGIN { GetDataAndResults }
Digits:=1-TRUNC(0.434294481*LN(Error));
REPEAT
GetReal('Enter distribution mean or [ENTER] to quit : ',Mean,Quit);
IF (NOT Quit) THEN
BEGIN
Getsd(sd,Quit);
IF (NOT Quit) THEN
BEGIN
GetReal('Enter lower limit : ',a,Quit);
IF (NOT Quit) THEN
BEGIN
GetReal('Enter upper limit : ',b,Quit);
IF (NOT Quit) THEN
BEGIN
WRITE('P[aóxób]=');
WRITELN(Probability((a-Mean)/sd,(b-Mean)/sd,Error):Digits+2:Digits)
END
END
END
END
UNTIL Quit
END; { GetDataAndResults }

BEGIN { Prob }
IF (ParamCount=0) THEN
GetDataAndResults(1.0E-6)
ELSE
IF (ParamCount>1) THEN
WRITELN('Only zero or one parameter should be used')
ELSE
BEGIN
ErrorString:=ParamStr(ParamCount);
Val(ErrorString,Error,Code);
IF (Code=0) THEN
BEGIN
IF (Error>0.1) THEN
WRITELN('Tolerance is too large')
ELSE
BEGIN
IF (Error<0.0) THEN
WRITELN('Tolerance must be positive')
ELSE
BEGIN
IF (Error<1.0E-16) THEN
WRITELN('Tolerance is too small')
ELSE
GetDataAndResults(Error)
END
END
END
ELSE
WRITELN('Invalid tolerance')

END
END. { Prob }