← Back to branch summary

~ubuntu-branches/ubuntu/utopic/mricron/utopic

« back to all changes in this revision

Viewing changes to fpmath/demo/stat/histo.pas

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-07-29 22:07:43 UTC
  • Revision ID: james.westby@ubuntu.com-20100729220743-q621ts2zj806gu0n
Tags: upstream-0.20100725.1~dfsg.1
Import upstream version 0.20100725.1~dfsg.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
fpmath/demo/stat/histo.pas
 
1
{ ******************************************************************
 
2
  Statistical distribution and Histogram
 
3
  ****************************************************************** }
 
4
 
 
5
program histo;
 
6
 
 
7
uses
 
8
  tpmath, tpgraph;
 
9
 
 
10
const
 
11
  N     = 30;    { Number of values }
 
12
  Alpha = 0.05;  { Significance level }
 
13
 
 
14
{ Hemoglobin concentrations in men }
 
15
const HbM : array[1..N] of Float =
 
16
(141, 144, 146, 148, 149, 150, 150, 151, 153, 153,
 
17
 153, 154, 155, 156, 156, 160, 160, 160, 163, 164,
 
18
 164, 165, 166, 168, 168, 170, 172, 172, 176, 179);
 
19
 
 
20
var
 
21
  M, S : Float;  { Mean and standard deviation }
 
22
 
 
23
function PltFunc(X : Float) : Float;
 
24
{ ------------------------------------------------------------------
 
25
  Function to be plotted (density of normal distribution)
 
26
  ------------------------------------------------------------------ }
 
27
 
 
28
begin
 
29
  PltFunc := DNorm((X - M) / S) / S;
 
30
end;
 
31
 
 
32
procedure WriteResults(C       : PStatClassVector;
 
33
                       Ncls    : Integer;
 
34
                       Calc    : PVector;
 
35
                       Khi2, G : Float;
 
36
                       DoF     : Integer;
 
37
                       K2c     : Float);
 
38
{ ------------------------------------------------------------------
 
39
  Writes results to screen
 
40
  ------------------------------------------------------------------ }
 
41
 
 
42
const
 
43
  Line1 = '-----------------------------';
 
44
 
 
45
var
 
46
  Sum : Float;
 
47
  I   : Integer;
 
48
 
 
49
begin
 
50
  Writeln('Statistical distribution');
 
51
 
 
52
  Writeln(Line1);
 
53
  Writeln('     Inf    Sup      N  Ncalc');
 
54
  Writeln(Line1);
 
55
 
 
56
  Sum := 0.0;
 
57
  for I := 1 to Ncls do
 
58
    begin
 
59
      Writeln(C^[I].Inf:8:0, C^[I].Sup:7:0, C^[I].N:7, Calc^[I]:7:2);
 
60
      Sum := Sum + Calc^[I];
 
61
    end;
 
62
 
 
63
  Writeln(Line1);
 
64
  Writeln('Total            ', N:5, Sum:7:2);
 
65
  Writeln(Line1);
 
66
 
 
67
  Writeln;
 
68
  Writeln('Comparison with normal distribution:');
 
69
  Writeln;
 
70
 
 
71
  Writeln('Pearson''s Khi-2           = ', Khi2:10:4, '  (', DoF, ' DoF)');
 
72
  Writeln('Woolf''s G                 = ', G:10:4, '  (', DoF, ' DoF)');
 
73
  Writeln('Critical value (p = ', Alpha:4:2, ') = ', K2c:10:4);
 
74
end;
 
75
 
 
76
procedure PlotGraph(C                  : PStatClassVector;
 
77
                    Ncls               : Integer;
 
78
                    Xmin , Xmax, Xstep : Float);
 
79
{ ------------------------------------------------------------------
 
80
  Plots histogram and normal curve
 
81
  ------------------------------------------------------------------ }
 
82
 
 
83
var
 
84
  Ymin, Ymax, Ystep : Float;    { Oy scale }
 
85
  Npts              : Integer;  { Number of points }
 
86
  X, Y              : PVector;  { Point coordinates }
 
87
  I, J              : Integer;  { Loop variables }
 
88
 
 
89
begin
 
90
  if not InitGraphics(9, 2, 'c:\tp\bgi') then  { 640x480 16 color }
 
91
    begin
 
92
      Writeln('Unable to set graphic mode');
 
93
      Exit;
 
94
    end;
 
95
 
 
96
  SetWindow(15, 85, 15, 85, True);
 
97
 
 
98
  { The histogram is plotted as a continuous curve.
 
99
    Each bar of the histogram is defined by 4 points }
 
100
 
 
101
  Npts := 4 * Ncls;
 
102
 
 
103
  DimVector(X, Npts);
 
104
  DimVector(Y, Npts);
 
105
 
 
106
  for I := 1 to Ncls do
 
107
    with C^[I] do
 
108
      begin
 
109
        J := 4 * (I - 1) + 1; X^[J] := Inf;  { Y^[J] := 0 }
 
110
        Inc(J); X^[J] := Inf; Y^[J] := D;
 
111
        Inc(J); X^[J] := Sup; Y^[J] := D;
 
112
        Inc(J); X^[J] := Sup;                { Y^[J] := 0 }
 
113
      end;
 
114
 
 
115
  { Set scale on Oy, making sure that it starts from 0 }
 
116
 
 
117
  AutoScale(Y, 1, Ncls, LinScale, Ymin, Ymax, Ystep);
 
118
 
 
119
  if Ymin <> 0.0 then Ymin := 0.0;
 
120
  Ymax := Ymax + Ystep;
 
121
 
 
122
  SetOxScale(LinScale, Xmin, Xmax, Xstep);
 
123
  SetOyScale(LinScale, Ymin, Ymax, Ystep);
 
124
 
 
125
  SetGraphTitle('Statistical Distribution');
 
126
  SetOxTitle('X');
 
127
  SetOyTitle('Frequency Density');
 
128
 
 
129
  PlotOxAxis;
 
130
  PlotOyAxis;
 
131
 
 
132
  WriteGraphTitle;
 
133
 
 
134
  SetClipping(True);
 
135
 
 
136
  { Plot histogram and normal curve }
 
137
 
 
138
  SetPointParam(1, 0, 0, 0);  { Don't show points on histogram }
 
139
  SetLineParam(1, 1, 3, 1);   { Use thick lines }
 
140
 
 
141
  PlotCurve(X, Y, 1, Npts, 1);
 
142
 
 
143
  PlotFunc({$IFDEF FPC}@{$ENDIF}PltFunc, Xmin, Xmax, 2);
 
144
 
 
145
  DelVector(X, Npts);
 
146
  DelVector(Y, Npts);
 
147
 
 
148
  Readln;
 
149
 
 
150
  LeaveGraphics;
 
151
end;
 
152
 
 
153
{ ******************************************************************
 
154
  Main program
 
155
  ****************************************************************** }
 
156
 
 
157
var
 
158
  X     : PVector;           { Data }
 
159
  C     : PStatClassVector;  { Statistical classes }
 
160
  Ncls  : Integer;           { Number of classes }
 
161
  Obs   : PIntVector;        { Observed frequencies }
 
162
  Calc  : PVector;           { Calculated frequencies }
 
163
  Khi2  : Float;             { Pearson's Khi-2 }
 
164
  G     : Float;             { Woolf's G }
 
165
  K2c   : Float;             { Theoretical Khi-2 }
 
166
  DoF   : Integer;           { Degrees of freedom }
 
167
  T     : Float;             { Standard normal variable }
 
168
  F0, F : Float;             { Cumulative probability }
 
169
  XMin,
 
170
  XMax,
 
171
  XStep : Float;             { Scale on Ox }
 
172
  I     : Integer;           { Loop variable }
 
173
 
 
174
begin
 
175
  { Read data }
 
176
  DimVector(X, N);
 
177
  for I := 1 to N do
 
178
    X^[I] := HbM[I];
 
179
 
 
180
  { Sort data if necessary }
 
181
  { QSort(X, 1, N);        }
 
182
 
 
183
  { Compute an appropriate interval for the set of values }
 
184
  Interval(X^[1], X^[N], 5, 10, XMin, XMax, XStep);
 
185
 
 
186
  { Compute number of classes and dimension arrays }
 
187
  Ncls := Round((Xmax - Xmin) / XStep);
 
188
 
 
189
  DimStatClassVector(C, Ncls);
 
190
  DimIntVector(Obs, Ncls);
 
191
  DimVector(Calc, Ncls);
 
192
 
 
193
  { Compute distribution }
 
194
  Distrib(X, 1, N, Xmin, Xmax, XStep, C);
 
195
 
 
196
  { Compute mean and S.D. }
 
197
  M := Mean(X, 1, N);
 
198
  S := StDev(X, 1, N, M);
 
199
 
 
200
  { Compute theoretical values }
 
201
  F0 := 0.0;
 
202
  for I := 1 to Ncls do
 
203
    begin
 
204
      if I = Ncls then
 
205
        F := 1.0
 
206
      else
 
207
        begin
 
208
          T := (C^[I].Sup - M) / S;
 
209
          F := FNorm(T);
 
210
        end;
 
211
      Calc^[I] := N * (F - F0);
 
212
      Obs^[I] := C^[I].N;
 
213
      F0 := F;
 
214
    end;
 
215
 
 
216
  { Perform Khi-2 and Woolf tests }
 
217
  Khi2_Conform(Ncls, 2, Obs, Calc, Khi2, DoF);
 
218
  Woolf_Conform(Ncls, 2, Obs, Calc, G, DoF);
 
219
 
 
220
  { Compute critical value }
 
221
  K2c := InvKhi2(DoF, 1.0 - Alpha);
 
222
 
 
223
  { Print results }
 
224
  WriteResults(C, Ncls, Calc, Khi2, G, DoF, K2c);
 
225
 
 
226
  Readln;
 
227
 
 
228
  { Plot histogram }
 
229
  PlotGraph(C, Ncls, Xmin - Xstep, Xmax + Xstep, Xstep);
 
230
end.