← Back to branch summary

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

« back to all changes in this revision

Viewing changes to fpmath/demo/curfit/reglin.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/curfit/reglin.pas
 
1
{ ******************************************************************
 
2
  This program performs a least squares fit of a straight line:
 
3
                         Y = B(0) + B(1) * X
 
4
  ****************************************************************** }
 
5
 
 
6
program reglin;
 
7
 
 
8
uses
 
9
  tpmath, tpgraph;
 
10
 
 
11
const
 
12
  N     = 5;     { Number of points }
 
13
  Alpha = 0.05;  { Significance level }
 
14
 
 
15
{ Data }
 
16
const
 
17
  X : array[1..N] of Float = (10,     20,     30,    40,     50);
 
18
  Y : array[1..N] of Float = ( 0.1865, 0.3616, 0.537, 0.7359, 0.9238);
 
19
 
 
20
var
 
21
  B : PVector;  { Regression parameters }
 
22
 
 
23
function PltFunc(X : Float) : Float;
 
24
{ ------------------------------------------------------------------
 
25
  Function to be plotted
 
26
  ------------------------------------------------------------------ }
 
27
 
 
28
begin
 
29
  PltFunc := B^[0] + B^[1] * X
 
30
end;
 
31
 
 
32
procedure WriteResults(X, Y, Ycalc, B : PVector;
 
33
                       V              : PMatrix;
 
34
                       Test           : TRegTest;
 
35
                       Tc, Fc         : Float);
 
36
{ ------------------------------------------------------------------
 
37
  Writes results to screen
 
38
  ------------------------------------------------------------------ }
 
39
 
 
40
var
 
41
  Line1,
 
42
  Line2 : String;   { Separating lines }
 
43
  Delta : Float;    { Residual }
 
44
  Sr    : Float;    { Residual standard deviation }
 
45
  SB    : Float;    { Standard deviations of parameters }
 
46
  I     : Integer;  { Loop variable }
 
47
 
 
48
begin
 
49
  Line1 := StrChar(73, '-');
 
50
  Line2 := StrChar(73, '=');
 
51
 
 
52
  WriteLn(Line2);
 
53
  WriteLn('Linear regression: Y = B(0) + B(1) * X');
 
54
  WriteLn(Line1);
 
55
 
 
56
  WriteLn('Parameter    Est.value         Std.dev.        ',
 
57
          (100 * (1 - Alpha)):2:0, '% Confidence Interval');
 
58
 
 
59
  WriteLn(Line1);
 
60
 
 
61
  for I := 0 to 1 do
 
62
    begin
 
63
      SB := Sqrt(V^[I]^[I]);
 
64
      WriteLn('B(', I:1, ')', B^[I]:17:8, SB:17:8,
 
65
              (B^[I] - Tc * SB):17:8, ';', (B^[I] + Tc * SB):17:8);
 
66
    end;
 
67
 
 
68
  WriteLn(Line1);
 
69
 
 
70
  WriteLn('Number of observations            : n           = ', N:5);
 
71
 
 
72
  with Test do
 
73
    begin
 
74
      Sr := Sqrt(Vr);
 
75
      WriteLn('Residual error                    : s           = ', Sr:10:4);
 
76
      WriteLn('Coefficient of correlation        : r           = ', (Sgn(B^[1]) * Sqrt(R2)):10:4);
 
77
      WriteLn('Coefficient of determination      : r2          = ', R2:10:4);
 
78
      WriteLn('Adjusted coeff. of determination  : r2a         = ', R2a:10:4);
 
79
      WriteLn('Variance ratio (explained/resid.) : F(', Nu1:3, ', ', Nu2:3, ') = ', F:10:4);
 
80
      WriteLn('Critical variance ratio           : F(p = ', (1 - Alpha):4:2, ') = ', Fc:10:4);
 
81
    end;
 
82
 
 
83
  WriteLn(Line1);
 
84
  WriteLn('  i        Y obs.       Y calc.      Residual      Std.dev.      Std.res.');
 
85
  WriteLn(Line1);
 
86
 
 
87
  for I := 1 to N do
 
88
    begin
 
89
      Delta := Y^[I] - Ycalc^[I];
 
90
      WriteLn(I:3, Y^[I]:14:4, Ycalc^[I]:14:4, Delta:14:4, Sr:14:4, (Delta / Sr):14:4);
 
91
    end;
 
92
 
 
93
  WriteLn(Line2);
 
94
end;
 
95
 
 
96
procedure PlotGraph(X, Y, B : PVector);
 
97
{ ------------------------------------------------------------------
 
98
  Plots histogram and normal curve
 
99
  ------------------------------------------------------------------ }
 
100
 
 
101
var
 
102
  Xmin, Xmax, Xstep : Float;    { Ox scale }
 
103
  Ymin, Ymax, Ystep : Float;    { Oy scale }
 
104
 
 
105
begin
 
106
  if not InitGraphics(9, 2, 'c:\tp\bgi') then  { 640x480 16 color }
 
107
    begin
 
108
      Writeln('Unable to set graphic mode');
 
109
      Exit;
 
110
    end;
 
111
 
 
112
  SetWindow(15, 85, 15, 85, True);
 
113
 
 
114
  AutoScale(X, 1, N, LinScale, Xmin, Xmax, Xstep);
 
115
  AutoScale(Y, 1, N, LinScale, Ymin, Ymax, Ystep);
 
116
 
 
117
  SetOxScale(LinScale, Xmin, Xmax, Xstep);
 
118
  SetOyScale(LinScale, Ymin, Ymax, Ystep);
 
119
 
 
120
  SetGraphTitle('Linear Regression');
 
121
  SetOxTitle('X');
 
122
  SetOyTitle('Y');
 
123
 
 
124
  PlotOxAxis;
 
125
  PlotOyAxis;
 
126
 
 
127
  WriteGraphTitle;
 
128
 
 
129
  SetClipping(True);
 
130
 
 
131
  SetLineParam(1, 0, 0, 0);  { Don't connect points }
 
132
  PlotCurve(X, Y, 1, N, 1);
 
133
 
 
134
  PlotFunc({$IFDEF FPC}@{$ENDIF}PltFunc, Xmin, Xmax, 2);
 
135
 
 
136
  Readln;
 
137
 
 
138
  LeaveGraphics;
 
139
end;
 
140
 
 
141
{ ******************************************************************
 
142
  Main program
 
143
  ****************************************************************** }
 
144
 
 
145
var
 
146
  XX, YY : PVector;   { Data }
 
147
  Ycalc  : PVector;   { Computed Y values }
 
148
  V      : PMatrix;   { Variance-covariance matrix }
 
149
  Test   : TRegTest;  { Statistical tests }
 
150
  Tc     : Float;     { Critical t value }
 
151
  Fc     : Float;     { Critical F value }
 
152
  I      : Integer;   { Loop variable }
 
153
 
 
154
begin
 
155
  { Dimension arrays }
 
156
  DimVector(XX, N);
 
157
  DimVector(YY, N);
 
158
  DimVector(Ycalc, N);
 
159
  DimVector(B, 1);
 
160
  DimMatrix(V, 1, 1);
 
161
 
 
162
  { Read data }
 
163
  for I := 1 to N do
 
164
    begin
 
165
      XX^[I] := X[I];
 
166
      YY^[I] := Y[I];
 
167
    end;
 
168
 
 
169
  { Perform regression }
 
170
  LinFit(XX, YY, 1, N, B, V);
 
171
 
 
172
  { Compute predicted Y values }
 
173
  for I := 1 to N do
 
174
    Ycalc^[I] := B^[0] + B^[1] * XX^[I];
 
175
 
 
176
  { Update variance-covariance matrix and compute statistical tests }
 
177
  RegTest(YY, Ycalc, 1, N, V, 0, 1, Test);
 
178
 
 
179
  { Compute Student's t and Snedecor's F }
 
180
  Tc := InvStudent(N - 2, 1 - 0.5 * Alpha);
 
181
  Fc := InvSnedecor(1, N - 2, 1 - Alpha);
 
182
 
 
183
  { Write results }
 
184
  WriteResults(XX, YY, Ycalc, B, V, Test, Tc, Fc);
 
185
 
 
186
  { Plot curve }
 
187
  PlotGraph(XX, YY, B);
 
188
end.