96 lines
2.4 KiB
Text
96 lines
2.4 KiB
Text
(* Statistics/Normal distribution *)
|
|
|
|
block
|
|
|
|
(* Generates normally distributed random numbers with mean 0 and standard deviation 1 *)
|
|
unit random_normal: function: real;
|
|
const
|
|
pi = 3.14159265358979;
|
|
begin
|
|
result := cos(2.0 * pi * random) * sqrt(-2.0 * ln(random))
|
|
end random_normal;
|
|
|
|
unit normal_stats: procedure (sample_size: integer);
|
|
var
|
|
r: arrayof real,
|
|
sum, mean, diff, scale, sd: real,
|
|
h: arrayof integer,
|
|
h_sum, adj, i, j: integer;
|
|
begin
|
|
if sample_size < 1 then return fi;
|
|
|
|
array r dim (1 : sample_size);
|
|
array h dim (-1 : 10); (* all zero by default *)
|
|
sum := 0.0;
|
|
h_sum := 0;
|
|
|
|
(* Generate 'sample_size' normally distributed random numbers with mean 0.5
|
|
and standard deviation 0.25.
|
|
Calculate their sum and in which box they will fall when drawing
|
|
the histogram. *)
|
|
for i := 1 to sample_size
|
|
do
|
|
r(i) := .5 + random_normal / 4.0;
|
|
sum := sum + r(i);
|
|
if r(i) < 0.0
|
|
then h(-1) := h(-1) + 1
|
|
else
|
|
if r(i) >= 1.0
|
|
then h(10) := h(10) + 1
|
|
else h(entier(r(i) * 10)) := h(entier(r(i) * 10)) + 1
|
|
fi
|
|
fi
|
|
od;
|
|
|
|
for i := -1 to 10 do h_sum := h_sum + h(i) od;
|
|
(* Adjust one of the h() values if necessary to ensure h_sum = sample_size *)
|
|
adj := sample_size - h_sum;
|
|
if adj =/= 0
|
|
then
|
|
for i := -1 to 10
|
|
do
|
|
h(i) := h(i) + adj;
|
|
if h(i) >= 0 then exit fi;
|
|
h(i) := h(i) - adj
|
|
od
|
|
fi;
|
|
mean := sum / sample_size;
|
|
|
|
(* Now calculate their standard deviation *)
|
|
sum := 0.0;
|
|
for i := 1 to sample_size
|
|
do
|
|
diff := r(i) - mean;
|
|
sum := sum + diff * diff
|
|
od;
|
|
sd := sqrt(sum / sample_size);
|
|
|
|
(* Draw a histogram of the data with interval 0.1
|
|
If sample size > 300 then normalize histogram to 300 *)
|
|
scale := 1.0;
|
|
if sample_size > 300 then scale := 300.0 / sample_size fi;
|
|
writeln("Sample size ", sample_size);
|
|
writeln(" Mean ", mean: 8: 6, " SD ", sd: 8: 6);
|
|
for i := -1 to 10
|
|
do
|
|
if i = -1 then write("< 0.00 : ")
|
|
else
|
|
if i = 10
|
|
then write(">=1.00 : ")
|
|
else write(" ", i / 10.0 : 4: 2, " : ")
|
|
fi
|
|
fi;
|
|
write(h(i): 5, " ");
|
|
for j := 1 to entier(h(i) * scale + .5) do write("*") od;
|
|
writeln
|
|
od;
|
|
end normal_stats;
|
|
|
|
begin
|
|
call normal_stats(100);
|
|
writeln;
|
|
call normal_stats(1000);
|
|
writeln;
|
|
call normal_stats(10000);
|
|
writeln;
|
|
end;
|