1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
|
{ ******************************************************************
Multiple linear regression (Singular Value Decomposition)
****************************************************************** }
unit usvdfit;
interface
uses
utypes, usvd;
procedure SVDFit(X : PMatrix;
Y : PVector;
Lb, Ub, Nvar : Integer;
ConsTerm : Boolean;
SVDTol : Float;
B : PVector;
V : PMatrix);
{ ------------------------------------------------------------------
Multiple linear regression: Y = B(0) + B(1) * X + B(2) * X2 + ...
------------------------------------------------------------------
Input parameters: X = matrix of independent variables
Y = vector of dependent variable
Lb, Ub = array bounds
Nvar = number of independent variables
ConsTerm = presence of constant term B(0)
Output parameters: B = regression parameters
V = inverse matrix
------------------------------------------------------------------ }
procedure WSVDFit(X : PMatrix;
Y, S : PVector;
Lb, Ub, Nvar : Integer;
ConsTerm : Boolean;
SVDTol : Float;
B : PVector;
V : PMatrix);
{ ----------------------------------------------------------------------
Weighted multiple linear regression
----------------------------------------------------------------------
S = standard deviations of observations
Other parameters as in SVDFit
---------------------------------------------------------------------- }
implementation
type
TRegMode = (UNWEIGHTED, WEIGHTED);
procedure GenSVDFit(Mode : TRegMode;
X : PMatrix;
Y, S : PVector;
Lb, Ub, Nvar : Integer;
ConsTerm : Boolean;
SVDTol : Float;
B : PVector;
V : PMatrix);
{ ----------------------------------------------------------------------
General multiple linear regression routine (SVD algorithm)
---------------------------------------------------------------------- }
var
U : PMatrix; { Matrix of independent variables for SVD }
Z : PVector; { Vector of dependent variables for SVD }
S1 : PVector; { Singular values }
S2inv : PVector; { Inverses of squared singular values }
V1 : PMatrix; { Orthogonal matrix from SVD }
LbU : Integer; { Lower bound of U matrix in both dim. }
UbU : Integer; { Upper bound of U matrix in 1st dim. }
I, J, K : Integer; { Loop variables }
Sigma : Float; { Square root of weight }
Sum : Float; { Element of variance-covariance matrix }
begin
if Ub - Lb < Nvar then
begin
SetErrCode(MatErrDim);
Exit;
end;
if Mode = WEIGHTED then
for K := Lb to Ub do
if S^[K] <= 0.0 then
begin
SetErrCode(MatSing);
Exit;
end;
{ ----------------------------------------------------------
Prepare arrays for SVD :
If constant term, use U[0..(N - Lb), 0..Nvar]
and Z[0..(N - Lb)]
else use U[1..(N - Lb + 1), 1..Nvar]
and Z[1..(N - Lb + 1)]
since the lower bounds of U for the SVD routine must be
the same in both dimensions
---------------------------------------------------------- }
if ConsTerm then
begin
LbU := 0;
UbU := Ub - Lb;
end
else
begin
LbU := 1;
UbU := Ub - Lb + 1;
end;
{ Dimension arrays }
DimMatrix(U, UbU, Nvar);
DimVector(Z, UbU);
DimVector(S1, Nvar);
DimVector(S2inv, Nvar);
DimMatrix(V1, Nvar, Nvar);
if Mode = UNWEIGHTED then
for I := LbU to UbU do
begin
K := I - LbU + Lb;
Z^[I] := Y^[K];
if ConsTerm then
U^[I]^[0] := 1.0;
for J := 1 to Nvar do
U^[I]^[J] := X^[K]^[J];
end
else
for I := LbU to UbU do
begin
K := I - LbU + Lb;
Sigma := 1.0 / S^[K];
Z^[I] := Y^[K] * Sigma;
if ConsTerm then
U^[I]^[0] := Sigma;
for J := 1 to Nvar do
U^[I]^[J] := X^[K]^[J] * Sigma;
end;
{ Perform singular value decomposition }
SV_Decomp(U, LbU, UbU, Nvar, S1, V1);
if MathErr = MatOk then
begin
{ Set the lowest singular values to zero }
SV_SetZero(S1, LbU, Nvar, SVDTol);
{ Solve the system }
SV_Solve(U, S1, V1, Z, LbU, UbU, Nvar, B);
{ Compute variance-covariance matrix }
for I := LbU to Nvar do
if S1^[I] > 0.0 then
S2inv^[I] := 1.0 / Sqr(S1^[I])
else
S2inv^[I] := 0.0;
for I := LbU to Nvar do
for J := LbU to I do
begin
Sum := 0.0;
for K := LbU to Nvar do
Sum := Sum + V1^[I]^[K] * V1^[J]^[K] * S2inv^[K];
V^[I]^[J] := Sum;
V^[J]^[I] := Sum;
end;
end;
DelMatrix(U, UbU, Nvar);
DelVector(Z, UbU);
DelVector(S1, Nvar);
DelVector(S2inv, Nvar);
DelMatrix(V1, Nvar, Nvar);
end;
procedure SVDFit(X : PMatrix;
Y : PVector;
Lb, Ub, Nvar : Integer;
ConsTerm : Boolean;
SVDTol : Float;
B : PVector;
V : PMatrix);
var
S : PVector;
begin
S := nil;
GenSVDFit(UNWEIGHTED, X, Y, S, Lb, Ub, Nvar, ConsTerm, SVDTol, B, V);
end;
procedure WSVDFit(X : PMatrix;
Y, S : PVector;
Lb, Ub, Nvar : Integer;
ConsTerm : Boolean;
SVDTol : Float;
B : PVector;
V : PMatrix);
begin
GenSVDFit(WEIGHTED, X, Y, S, Lb, Ub, Nvar, ConsTerm, SVDTol, B, V);
end;
end.
|
