← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/mricron/oneiric

« back to all changes in this revision

Viewing changes to fpmath/umulfit.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/umulfit.pas
 
1
{ ******************************************************************
 
2
  Multiple linear regression (Gauss-Jordan method)
 
3
  ****************************************************************** }
 
4
 
 
5
unit umulfit;
 
6
 
 
7
interface
 
8
 
 
9
uses
 
10
  utypes, ulineq;
 
11
 
 
12
 
 
13
procedure MulFit(X            : PMatrix;
 
14
                 Y            : PVector;
 
15
                 Lb, Ub, Nvar : Integer;
 
16
                 ConsTerm     : Boolean;
 
17
                 B            : PVector;
 
18
                 V            : PMatrix);
 
19
{ ------------------------------------------------------------------
 
20
  Multiple linear regression: Y = B(0) + B(1) * X + B(2) * X2 + ...
 
21
  ------------------------------------------------------------------
 
22
  Input parameters:  X        = matrix of independent variables
 
23
                     Y        = vector of dependent variable
 
24
                     Lb, Ub   = array bounds
 
25
                     Nvar     = number of independent variables
 
26
                     ConsTerm = presence of constant term B(0)
 
27
  Output parameters: B        = regression parameters
 
28
                     V        = inverse matrix
 
29
  ------------------------------------------------------------------ }
 
30
 
 
31
procedure WMulFit(X            : PMatrix;
 
32
                  Y, S         : PVector;
 
33
                  Lb, Ub, Nvar : Integer;
 
34
                  ConsTerm     : Boolean;
 
35
                  B            : PVector;
 
36
                  V            : PMatrix);
 
37
{ ----------------------------------------------------------------------
 
38
  Weighted multiple linear regression
 
39
  ----------------------------------------------------------------------
 
40
  S = standard deviations of observations
 
41
  Other parameters as in MulFit
 
42
  ---------------------------------------------------------------------- }
 
43
 
 
44
implementation
 
45
 
 
46
procedure MulFit(X            : PMatrix;
 
47
                 Y            : PVector;
 
48
                 Lb, Ub, Nvar : Integer;
 
49
                 ConsTerm     : Boolean;
 
50
                 B            : PVector;
 
51
                 V            : PMatrix);
 
52
 
 
53
  var
 
54
    Lb1     : Integer;  { Index of first param. (0 if cst term, 1 otherwise) }
 
55
    I, J, K : Integer;  { Loop variables }
 
56
    Det     : Float;    { Determinant }
 
57
 
 
58
  begin
 
59
    if Ub - Lb < Nvar then
 
60
      begin
 
61
        SetErrCode(MatErrDim);
 
62
        Exit;
 
63
      end;
 
64
 
 
65
    { Initialize }
 
66
    for I := 0 to Nvar do
 
67
      begin
 
68
        for J := 0 to Nvar do
 
69
          V^[I]^[J] := 0.0;
 
70
        B^[I] := 0.0;
 
71
      end;
 
72
 
 
73
    { If constant term, set line 0 and column 0 of matrix V }
 
74
    if ConsTerm then
 
75
      begin
 
76
        V^[0]^[0] := Int(Ub - Lb + 1);
 
77
        for K := Lb to Ub do
 
78
          begin
 
79
            for J := 1 to Nvar do
 
80
              V^[0]^[J] := V^[0]^[J] + X^[K]^[J];
 
81
            B^[0] := B^[0] + Y^[K];
 
82
          end;
 
83
        for J := 1 to Nvar do
 
84
          V^[J]^[0] := V^[0]^[J];
 
85
      end;
 
86
 
 
87
    { Set other elements of V }
 
88
    for K := Lb to Ub do
 
89
      for I := 1 to Nvar do
 
90
        begin
 
91
          for J := I to Nvar do
 
92
            V^[I]^[J] := V^[I]^[J] + X^[K]^[I] * X^[K]^[J];
 
93
          B^[I] := B^[I] + X^[K]^[I] * Y^[K];
 
94
        end;
 
95
 
 
96
    { Fill in symmetric matrix }
 
97
    for I := 2 to Nvar do
 
98
      for J := 1 to Pred(I) do
 
99
        V^[I]^[J] := V^[J]^[I];
 
100
 
 
101
    { Solve normal equations }
 
102
    if ConsTerm then Lb1 := 0 else Lb1 := 1;
 
103
    LinEq(V, B, Lb1, Nvar, Det);
 
104
  end;
 
105
 
 
106
procedure WMulFit(X            : PMatrix;
 
107
                  Y, S         : PVector;
 
108
                  Lb, Ub, Nvar : Integer;
 
109
                  ConsTerm     : Boolean;
 
110
                  B            : PVector;
 
111
                  V            : PMatrix);
 
112
 
 
113
  var
 
114
    Lb1     : Integer;  { Index of first param. (0 if cst term, 1 otherwise) }
 
115
    I, J, K : Integer;  { Loop variables }
 
116
    W       : PVector;  { Vector of weights }
 
117
    WX      : Float;    { W * X }
 
118
    Det     : Float;    { Determinant }
 
119
 
 
120
  begin
 
121
    if Ub - Lb < Nvar then
 
122
      begin
 
123
        SetErrCode(MatErrDim);
 
124
        Exit;
 
125
      end;
 
126
 
 
127
    for K := Lb to Ub do
 
128
      if S^[K] <= 0.0 then
 
129
        begin
 
130
          SetErrCode(MatSing);
 
131
          Exit;
 
132
        end;
 
133
 
 
134
    DimVector(W, Ub);
 
135
 
 
136
    for K := Lb to Ub do
 
137
      W^[K] := 1.0 / Sqr(S^[K]);
 
138
 
 
139
    { Initialize }
 
140
    for I := 0 to Nvar do
 
141
      begin
 
142
        for J := 0 to Nvar do
 
143
          V^[I]^[J] := 0.0;
 
144
        B^[I] := 0.0;
 
145
      end;
 
146
 
 
147
    { If constant term, set line 0 and column 0 of matrix V }
 
148
    if ConsTerm then
 
149
      begin
 
150
        for K := Lb to Ub do
 
151
          begin
 
152
            V^[0]^[0] := V^[0]^[0] + W^[K];
 
153
            for J := 1 to Nvar do
 
154
              V^[0]^[J] := V^[0]^[J] + W^[K] * X^[K]^[J];
 
155
            B^[0] := B^[0] + W^[K] * Y^[K];
 
156
          end;
 
157
        for J := 1 to Nvar do
 
158
          V^[J]^[0] := V^[0]^[J];
 
159
      end;
 
160
 
 
161
    { Set other elements of V }
 
162
    for K := Lb to Ub do
 
163
      for I := 1 to Nvar do
 
164
        begin
 
165
          WX := W^[K] * X^[K]^[I];
 
166
          for J := I to Nvar do
 
167
            V^[I]^[J] := V^[I]^[J] + WX * X^[K]^[J];
 
168
          B^[I] := B^[I] + WX * Y^[K];
 
169
        end;
 
170
 
 
171
    { Fill in symmetric matrix }
 
172
    for I := 2 to Nvar do
 
173
      for J := 1 to Pred(I) do
 
174
        V^[I]^[J] := V^[J]^[I];
 
175
 
 
176
    { Solve normal equations }
 
177
    if ConsTerm then Lb1 := 0 else Lb1 := 1;
 
178
    LinEq(V, B, Lb1, Nvar, Det);
 
179
 
 
180
    DelVector(W, Ub);
 
181
  end;
 
182
 
 
183
end.
 
 
b'\\ No newline at end of file'