RosettaCodeData/Task/Statistics-Normal-distribution/Pascal/statistics-normal-distribut...

133 lines
3.3 KiB
ObjectPascal

Program Example40;
{$IFDEF FPC}
{$MOde objFPC}
{$ENDIF}
{ Program to demonstrate the randg function. }
Uses Math;
type
tTestData = extended;//because of math.randg
ttstfunc = function (mean, sd: tTestData): tTestData;
tExArray = Array of tTestData;
tSolution = record
SolExArr : tExArray;
SollowVal,
SolHighVal,
SolMean,
SolStdDiv : tTestData;
SolSmpCnt : LongInt;
end;
function getSol(genFunc:ttstfunc;Mean,StdDiv: tTestData;smpCnt: LongInt): tSolution;
var
GenValue,
sumValue,
sumsqrVal : extended;
Begin
with result do
Begin
SolSmpCnt := smpCnt;
SolMean := 0;
SolStdDiv := 0;
SolLowVal := Mean+50* StdDiv;
SolHighVal := Mean-50* StdDiv;
setlength(SolExArr,smpCnt);
if smpCnt <= 0 then
EXIT;
sumValue := 0;
sumsqrVal := 0;
repeat
GenValue := genFunc(Mean,StdDiv);
sumValue := sumvalue+GenValue;
sumsqrVal := sumsqrVal+sqr(GenValue);
IF GenValue < SollowVal then
SollowVal:= GenValue
else
IF GenValue > SolHighVal then
SolHighVal := GenValue;
dec(smpCnt);
SolExArr[smpCnt] := GenValue;
until smpCnt<= 0;
SolMean := sumValue/SolSmpCnt;
SolStdDiv := sqrt(sumsqrVal/SolSmpCnt-sqr(SolMean));
end;
end;
//http://wiki.freepascal.org/Generating_Random_Numbers#Normal_.28Gaussian.29_Distribution
function rnorm (mean, sd: tTestData): tTestData;
{Calculates Gaussian random numbers according to the Box-Müller approach}
var
u1, u2: extended;
begin
u1 := random;
u2 := random;
rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
end;
procedure Histo(const sol:TSolution;Colcnt,ColLen :LongInt);
var
CntHisto : array of integer;
LoLmt,HiLmt,span : tTestData;
i, j,cnt,maxCnt: LongInt;
sCross : Ansistring;
Begin
setlength(CntHisto,Colcnt);
with Sol do
Begin
span := solHighVal-solLowVal;
LoLmt := solLowVal;
writeln('Count: ',SolSmpCnt:10,' Mean ',SolMean:10:6,' StdDiv ',SolStdDIv:10:6);
writeln('span : ',span:10:5,' Low ',solLowVal:10:6,' high ',solHighVal:10:6);
end;
maxCnt := 0;
For j := 0 to Colcnt-1 do
Begin
HiLmt:= LoLmt+span/Colcnt;
cnt := 0;
with sol do
For i := 0 to High(SolExArr) do
IF (HiLmt > SolExArr[i]) AND (SolExArr[i]>= LoLmt) then
inc(cnt);
CntHisto[j] := cnt;
IF maxCnt < cnt then
maxCnt := cnt;
LoLmt:= HiLmt;
end;
inc(CntHisto[Colcnt]); // for HiLmt itself
writeln;
LoLmt := sol.solLowVal;
For i := 0 to Colcnt-1 do
Begin
Writeln(LoLmt:8:4,': ');
cnt:= Round(CntHisto[i]*ColLen/maxCnt);
setlength(sCross,cnt+3);
fillChar(sCross[1],3,' ');
fillChar(sCross[4],cnt,'#');
writeln(CntHisto[i]:10,sCross);
LoLmt := LoLmt+span/Colcnt;
end;
Writeln(sol.solHighVal:8:4,': ');
end;
const
cHistCnt = 11;
cColLen = 65;
cStdDiv = 0.25;
cMean = 20*cStdDiv;
var
mySol : tSolution;
begin
Randomize;
// test of randg of unit math
Writeln('function randg');
mySol := getSol(@randg,cMean,cMean*cStdDiv,100000);
Histo(mySol,cHistCnt,cColLen);
writeln;
// test of rnorm from wiki
Writeln('function rnorm');
mySol := getSol(@rnorm,cMean,cStdDiv,1000000);
Histo(mySol,cHistCnt,cColLen);
end.