pascal/lazarus-ccr
1
0
Fork 0
Code Issues Releases Activity
Files
e1c5977e0d1d9282ee080ff58bbda9f7dfab3d07
lazarus-ccr/applications/lazstats/source_orig/matrixlib.pas
wp_xxyyzz e1c5977e0d LazStats: Adding original source, part 6. 
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7885 8e941d3f-bd1b-0410-a28a-d453659cc2b4
2020-11-16 11:16:49 +00:00

2386 lines
79 KiB
ObjectPascal
Raw Blame History

unit MatrixLib;
{$mode objfpc}{$H+}
interface
uses
Classes, SysUtils, Globals, DictionaryUnit, OutPutUnit, Dialogs,
Math, FunctionsLib, DataProcs, MainUnit;
procedure GridVecProd(col1, col2 : integer; VAR Product : double; VAR Ngood : integer);
procedure GridXProd(NoSelected : integer;
VAR Selected : IntDyneVec;
VAR Product : DblDyneMat;
Augment : boolean;
VAR Ngood : integer);
procedure GridCovar(NoSelected : integer;
VAR Selected : IntDyneVec;
VAR Covar : DblDyneMat;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR errorcode : boolean;
VAR Ngood : integer);
procedure Correlations(NoSelected : integer;
VAR Selected : IntDyneVec;
VAR Correlations : DblDyneMat;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR errorcode : boolean;
VAR Ngood : integer);
procedure MATAxB(VAR a : DblDyneMat;
VAR b : DblDyneMat;
VAR c : DblDyneMat;
brows, bcols, crows, ccols : integer;
VAR errorcode : boolean);
function MATTRN(VAR a : DblDyneMat;
VAR b : DblDyneMat;
brows, bcols : integer): boolean;
procedure nonsymroots(a : DblDyneMat; nv : integer;
var nf : integer; c : real;
var v : DblDyneMat; var e : DblDyneVec;
var px : DblDyneVec;
var t : double;
var ev : double);
PROCEDURE ludcmp(VAR a: DblDyneMat; n: integer; VAR indx: IntDyneVec; VAR d: double);
procedure DETERM(VAR a : DblDyneMat; rows, cols : integer; VAR determ : double;
VAR errorcode : boolean);
procedure EffectCode(GridCol, min, max : integer;
FactLetter : string;
VAR startcol : integer;
VAR endcol : integer;
VAR novectors : integer);
procedure MReg(NoIndep : integer;
VAR IndepCols : IntDyneVec;
DepCol : integer;
VAR RowLabels : StrDyneVec;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR BWeights : DblDyneVec;
VAR BetaWeights : DblDyneVec;
VAR BStdErrs : DblDyneVec;
VAR Bttests : DblDyneVec;
VAR tProbs : DblDyneVec;
VAR R2 : double;
VAR stderrest : double;
VAR NCases : integer;
VAR errorcode : boolean;
PrintAll : boolean);
procedure Dynnonsymroots(var a : DblDyneMat; nv : integer;
var nf : integer; c : real;
var v : DblDyneMat; var e : DblDyneVec;
var px : DblDyneVec;
var t : double;
var ev : double);
function DynCorrelations(novars : integer;
VAR ColSelected : IntDyneVec;
VAR DataGrid : DblDyneMat;
VAR rmatrix : DblDyneMat;
VAR means : DblDyneVec;
VAR vars : DblDyneVec;
VAR stddevs : DblDyneVec;
NCases : integer;
ReturnType : integer) : integer;
procedure Predict(VAR ColNoSelected : IntDyneVec;
NoVars : integer;
VAR IndepInverse : DblDyneMat;
VAR Means : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR BetaWeights : DblDyneVec;
StdErrEst : double;
VAR IndepIndex : IntDyneVec;
NoIndepVars : integer);
procedure MReg2(NCases : integer;
NoVars : integer;
VAR NoIndepVars : integer;
VAR IndepIndex : IntDyneVec;
VAR corrs : DblDyneMat;
VAR IndepCorrs : DblDyneMat;
VAR RowLabels : StrDyneVec;
VAR R2 : double;
VAR BetaWeights : DblDyneVec;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR errorcode : integer;
VAR StdErrEst : double;
VAR constant : double;
probout : double;
Printit : boolean;
TestOut : boolean;
PrintInv : boolean);
procedure MATSUB(VAR a, b, c : DblDyneMat;
brows, bcols, crows, ccols : integer; VAR errorcode : boolean);
procedure IntArrayPrint(mat : IntDyneMat;
rows, cols : integer;
ytitle : string;
RowLabels, ColLabels : StrDyneVec;
Title : string);
procedure eigens(VAR a: DblDyneMat; Var d : DblDyneVec; n : integer);
PROCEDURE tred2(VAR a: DblDyneMat; n: integer; VAR d,e: DblDyneVec);
PROCEDURE tqli(VAR d,e: DblDyneVec; n: integer; VAR z: DblDyneMat);
function SEVS(nv,nf : integer;
c : double;
var r : DblDyneMat;
VAR v : DblDyneMat;
VAR e : DblDyneVec;
var p : DblDyneVec;
VAR nd : integer) : integer ;
function SCPF(VAR x,y : DblDyneMat; kx,ky,n,nd : integer) : double;
procedure MAT_PRINT(VAR xmat : DblDyneMat ;
rows, cols : integer;
VAR title : string;
VAR RowLabels: StrDyneVec;
VAR ColLabels: StrDyneVec;
NCases: integer);
procedure DynVectorPrint(VAR avector: DblDyneVec;
NoVars : integer;
title : string;
VAR Labels : StrDyneVec;
NCases : integer);
procedure scatplot(var x : DblDyneVec;
var y : DblDyneVec;
nocases : integer;
titlestr : string;
x_axis, y_axis : string;
x_min, x_max, y_min, y_max : double;
VAR VarLabels : StrDyneVec);
procedure DynIntMatPrint(mat : IntDyneMat;
rows, cols : integer;
ytitle : string;
RowLabels, ColLabels : StrDyneVec;
Title : string);
procedure SymMatRoots(A : DblDyneMat; M : integer; VAR E : DblDyneVec; VAR V : DblDyneMat);
procedure matinv(a, vtimesw, v, w: DblDyneMat; n: integer);
implementation
procedure GridVecProd(col1, col2 : integer; VAR Product : double; VAR Ngood : integer);
// Get the cross-product of two vectors
// col1 and col2 are grid columns
// Product is the vector product
// Ngood are the number of elements in the product not missing or filtered
var
i : integer;
Selected : IntDyneVec;
X1, X2 : double;
begin
SetLength(Selected,2);
Product := 0.0;
Selected[0] := col1;
Selected[1] := col2;
for i := 1 to NoCases do
begin
if NOT GoodRecord(i,2,Selected) then continue;
Ngood := Ngood + 1;
X1 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[col1,i]));
X2 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[col2,i]));
Product := Product + (X1 * X2);
end;
Selected := nil;
end;
//-------------------------------------------------------------------
procedure GridXProd(NoSelected : integer;
VAR Selected : IntDyneVec;
VAR Product : DblDyneMat;
Augment : boolean;
VAR Ngood : integer);
// Matrix product of a grid matrix and its transpose
// Product contains the cross-products matrix upon return
// Selected is a integer vector of grid columns of the vectors
// NoSelected is an integer of the number of grid vectors selected
// Ngood is the number of elements in a vector product not missing or filtered
// Augment is true if the augment matrix is to be obtained and is required
// to obtain means, variances, standard deviations in the correlation procedure
// and GridCovar procedure
var
i, j, k : integer;
Col1, Col2 : integer;
X1 : double;
Prod : double;
NoVars : integer;
N : double;
begin
// initialize
N := 0.0;
NoVars := 0;
for i := 1 to NoSelected do
for j := 1 to NoSelected do
Product[i-1,j-1] := 0.0;
if Augment then
begin
NoVars := NoSelected + 1;
for i := 1 to NoVars do
begin
Product[i-1,NoVars-1] := 0.0;
Product[NoVars-1,i-1] := 0.0;
end;
end;
// Do cross-products without augmentation
for i := 1 to NoSelected do // pre-matrix row (Grid transpose)
begin
for j := 1 to NoSelected do // post-matrix column (Grid)
begin
Ngood := 0;
Col1 := Selected[i-1];
Col2 := Selected[j-1];
GridVecProd(Col1,Col2,Prod, Ngood);
Product[i-1,j-1] := Prod;
end;
end;
if Augment then // do last column and row for augmented matrix
begin
for j := 1 to NoSelected do
begin
Col1 := Selected[j-1];
for k := 1 to NoCases do
begin
if NOT GoodRecord(k,NoSelected,Selected) then continue;
X1 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[Col1,k]));
Product[NoVars-1,j-1] := Product[NoVars-1,j-1] + X1;
Product[j-1,NoVars-1] := Product[j-1,NoVars-1] + X1;
end;
end;
for i := 1 to NoCases do // last cell of augmented matix
begin
if NOT GoodRecord(i,NoSelected,Selected) then continue;
N := N + 1.0;
end;
Product[NoVars-1,NoVars-1] := N;
Ngood := round(N);
end;
end;
//-------------------------------------------------------------------
procedure GridCovar(NoSelected : integer;
VAR Selected : IntDyneVec;
VAR Covar : DblDyneMat;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR errorcode : boolean;
VAR Ngood : integer);
// Obtains the variance/covariance matrix of variables in the grid
// NoSelected is the number of variables selected from the grid
// Selected is a vector of integers for the grid columns of selected variables
// Covar is the variance/covariance matrix returned
// Means, StdDevs, Variances are double vectors obtained from the augmented matrix
// errorcode is true if an error occurs due to 0 variance
// Ngood is the number of records in the cross-product of vectors
// This procedure calls the GridXProd procedure with augmentation true
// in order to obtain the means, variances and standard deviations
var
i, j : integer;
N : double;
Augment : boolean;
begin
// initialize
errorcode := false;
for i := 1 to NoSelected do
begin
Means[i-1] := 0.0;
Variances[i-1] := 0.0;
StdDevs[i-1] := 0.0;
end;
Augment := true; // augment to get intercept, means, variances, std.devs.
// get cross-products
GridXProd(NoSelected,Selected,Covar,Augment,Ngood);
// Get no. of records in cross-products
N := Ngood;
// Sums of squares are in diagonal, cross-products in off-diagonal cells
// Sums of X's are in the augmented column
// Get means and standard deviations first
for i := 1 to NoSelected do
begin
Means[i-1] := Covar[i-1,NoSelected] / N;
Variances[i-1] := Covar[i-1,i-1] - (Sqr(Covar[i-1,NoSelected]) / N);
Variances[i-1] := Variances[i-1] / (N - 1.0);
if Variances[i-1] > 0.0 then
StdDevs[i-1] := sqrt(Variances[i-1])
else
begin
StdDevs[i-1] := 0.0;
errorcode := true;
end;
end;
// Now get covariances
for i := 1 to NoSelected do
begin
for j := 1 to NoSelected do
begin
Covar[i-1,j-1] := Covar[i-1,j-1] -
((Covar[i-1,NoSelected] * Covar[j-1,NoSelected]) / N);
Covar[i-1,j-1] := Covar[i-1,j-1] / (N - 1);
end;
end;
end;
//-------------------------------------------------------------------
procedure Correlations(NoSelected : integer;
VAR Selected : IntDyneVec;
VAR Correlations : DblDyneMat;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR errorcode : boolean;
VAR Ngood : integer);
// Obtains the correlation matrix among grid variables
// NoSelected is the no. of grid variables selected for analysis
// Selected is a vector of integers of the grid variable columns selected
// Correlations are returned in the Correlations matrix
// Means, Variances, StdDevs are returned as double vectors
// errorcode is true if a 0 variance is detected
// Ngood is the number cases that do not contain missing values or are filtered
// This procedure calls the GridCovar procedure
var
i, j : integer;
begin
// get covariance matrix, means and standard deviations
GridCovar(NoSelected,Selected,Correlations,Means,Variances,StdDevs,errorcode, Ngood);
for i := 1 to NoSelected do
begin
for j := 1 to NoSelected do
begin
if (StdDevs[i-1] > 0.0) and (StdDevs[j-1] > 0.0) then
Correlations[i-1,j-1] := Correlations[i-1,j-1] /
(StdDevs[i-1] * StdDevs[j-1])
else
begin
Correlations[i-1,j-1] := 0.0;
errorcode := true;
end;
end;
end;
end;
//-------------------------------------------------------------------
procedure MATAxB(VAR a : DblDyneMat;
VAR b : DblDyneMat;
VAR c : DblDyneMat;
brows, bcols, crows, ccols : integer;
VAR errorcode : boolean);
// Product of matrix b times c with results returned in a
var i, j, k : integer;
begin
errorcode := FALSE;
if (bcols <> crows) then errorcode := TRUE
else
begin
for i := 0 to brows-1 do
begin
for j := 0 to ccols-1 do
begin
a[i,j] := 0.0;
for k := 0 to crows-1 do
a[i,j] := a[i,j] + b[i,k] * c[k,j];
end;
end;
end;
end; { of MATAxB }
//-------------------------------------------------------------------
function MATTRN(VAR a : DblDyneMat;
VAR b : DblDyneMat;
brows, bcols : integer): boolean;
// transpose the b matrix and return it in a
var i, j : integer;
errorcode : boolean;
begin
errorcode := FALSE;
if ((brows = 1) and (bcols = 1)) then errorcode := TRUE else
begin
for i := 0 to brows-1 do
for j := 0 to bcols-1 do
a[j,i] := b[i,j];
end;
MATTRN := errorcode;
end; { of mattrn }
//-------------------------------------------------------------------
procedure nonsymroots(a : DblDyneMat; nv : integer;
var nf : integer; c : real;
var v : DblDyneMat; var e : DblDyneVec;
var px : DblDyneVec;
var t : double;
var ev : double);
{ roots and vectors of a non symetric matrix. a is square matrix entered
and is destroyed in process. nv is number of variables (rows and columns )
of a. nf is the number of factorsto be extracted - is output as the number
which exceeded c, the minimum eigenvalue to be extracted. v is the output
matrix of column vectors of loadings. e is the output vector of roots. px
is the percentages of trace for factors. t is the trace of the matrix and
ev is the percent of trace extracted }
label 40;
var
y, z : DblDyneVec;
ek, e2, d : real;
i, j, k, m : integer;
begin
SetLength(y,nv);
SetLength(z,nv);
t := 0.0;
for i := 0 to nv-1 do t := t + a[i,i];
for k := 0 to nf-1 do
begin
for i := 0 to nv-1 do
begin
px[i] := 1.0;
y[i] := 1.0;
end;
e[k] := 1.0;
ek := 1.0;
for m := 1 to 25 do
begin
for i := 0 to nv-1 do
begin
v[i,k] := px[i] / e[k];
z[i] := y[i] / ek;
end;
for i := 0 to nv-1 do
begin
px[i] := 0.0;
for j := 0 to nv-1 do px[i] := px[i] + a[i,j] * v[j,k];
y[i] := 0.0;
for j := 0 to nv-1 do y[i] := y[i] + a[j,i] * z[j];
end;
e2 := 0.0;
for j := 0 to nv-1 do e2 := e2 + px[j] * v[j,k];
e[k] := sqrt(abs(e2));
ek := 0.0;
for j := 0 to nv-1 do ek := ek + y[j] * z[j];
ek := sqrt(abs(ek));
end;
if e2 >= sqr(c) then
begin
d := 0.0;
for j := 0 to nv-1 do d := d + v[j,k] * z[j];
d := e[k] / d;
for i := 0 to nv-1 do
for j := 0 to nv-1 do
a[i,j] := a[i,j] - v[i,k] * z[j] * d;
end
else begin
nf := k - 1;
goto 40;
end;
end;
40 : for i := 0 to nf-1 do px[i] := e[i] / t * 100.0;
ev := 0.0;
for i := 0 to nf-1 do ev := ev + px[i];
z := nil;
y := nil;
end; { of procedure nonsymroots }
//-------------------------------------------------------------------
PROCEDURE ludcmp(VAR a: DblDyneMat; n: integer; VAR indx: IntDyneVec; VAR d: double);
CONST tiny=1.0e-20;
VAR k,j,imax,i: integer;
sum,dum,big: double;
vv: DblDyneVec;
BEGIN
SetLength(vv,n);
d := 1.0; imax := 0;
FOR i := 1 to n DO BEGIN
big := 0.0;
FOR j := 1 to n DO IF (abs(a[i-1,j-1]) > big) THEN big := abs(a[i-1,j-1]);
IF (big = 0.0) THEN BEGIN
ShowMessage('Singular matrix in Lower-Upper Decomposition routine');
exit;
END;
vv[i-1] := 1.0/big
END;
FOR j := 1 to n DO BEGIN
IF (j > 1) THEN BEGIN
FOR i := 1 to j-1 DO BEGIN
sum := a[i-1,j-1];
IF (i > 1) THEN BEGIN
FOR k := 1 to i-1 DO BEGIN
sum := sum - a[i-1,k-1] * a[k-1,j-1]
END;
a[i-1,j-1] := sum
END
END
END;
big := 0.0;
FOR i := j to n DO BEGIN
sum := a[i-1,j-1];
IF (j > 1) THEN BEGIN
FOR k := 1 to j-1 DO BEGIN
sum := sum - a[i-1,k-1] * a[k-1,j-1]
END;
a[i-1,j-1] := sum
END;
dum := vv[i-1] * abs(sum);
IF (dum > big) THEN BEGIN
big := dum;
imax := i
END
END;
IF (j <> imax) THEN BEGIN
FOR k := 1 to n DO BEGIN
dum := a[imax-1,k-1];
a[imax-1,k-1] := a[j-1,k-1];
a[j-1,k-1] := dum
END;
d := -d;
vv[imax-1] := vv[j-1]
END;
indx[j-1] := imax;
IF (j <> n) THEN BEGIN
IF (a[j-1,j-1] = 0.0) THEN a[j-1,j-1] := tiny;
dum := 1.0/a[j-1,j-1];
FOR i := j+1 to n DO BEGIN
a[i-1,j-1] := a[i-1,j-1] * dum
END
END
END;
IF (a[n-1,n-1] = 0.0) THEN a[n-1,n-1] := tiny;
vv := nil;
END;
//-------------------------------------------------------------------
procedure DETERM(VAR a : DblDyneMat; rows, cols : integer; VAR determ : double;
VAR errorcode : boolean);
var indx : IntDyneVec;
i : integer;
begin
SetLength(indx,rows);
errorcode := FALSE;
if (rows <> cols) then errorcode := TRUE else
begin
LUDCMP(a, rows, indx, determ);
for i := 1 to rows do
determ := determ * a[i-1,i-1];
end;
end; { of determ }
//-------------------------------------------------------------------
procedure EffectCode(GridCol, min, max : integer;
FactLetter : string;
VAR startcol : integer;
VAR endcol : integer;
VAR novectors : integer);
var
levels, i, j, grp, col, cval : integer;
coef : IntDyneMat;
labelstr : string;
begin
// Routine for creating coded vectors representing group membership
// for purposes of multiple regression effects of group membership
levels := max - min + 1;
SetLength(coef,levels,levels);
novectors := levels - 1;
startcol := NoVariables + 1;
endcol := startcol + novectors - 1;
// setup grid for additional columns
for i := 1 to levels - 1 do
begin
labelstr := FactLetter + IntToStr(i);
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := labelstr;
OS3MainFrm.DataGrid.Cells[col,0] := labelstr;
end;
// get coefficients for effect coding
for i := 1 to levels do // group code
begin
for j := 1 to levels - 1 do // vector code
begin
if i = j then coef[i-1,j-1] := 1;
if i = levels then coef[i-1,j-1] := -1;
if (i <> j) and (i <> levels) then coef[i-1,j-1] := 0;
end;
end;
// code the cases using coefficients above
col := NoVariables - (levels - 1);
for i := 1 to levels - 1 do
begin
col := col + 1;
for j := 1 to NoCases do
begin
grp := round(StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[GridCol,j])))
- min + 1; // subject group code
cval := coef[grp-1,i-1]; // vector code
OS3MainFrm.DataGrid.Cells[col,j] := IntToStr(cval);
end;
end;
// NoVariables := NoVariables + novectors;
OS3MainFrm.NoVarsEdit.Text := IntToStr(NoVariables);
coef := nil;
end;
//-------------------------------------------------------------------
procedure MReg(NoIndep : integer;
VAR IndepCols : IntDyneVec;
DepCol : integer;
VAR RowLabels : StrDyneVec;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR BWeights : DblDyneVec;
VAR BetaWeights : DblDyneVec;
VAR BStdErrs : DblDyneVec;
VAR Bttests : DblDyneVec;
VAR tProbs : DblDyneVec;
VAR R2 : double;
VAR stderrest : double;
VAR NCases : integer;
VAR errorcode : boolean;
PrintAll : boolean);
var
i, j, k, N : integer;
X : DblDyneMat;
XT : DblDyneMat;
XTX : DblDyneMat;
XTY : DblDyneVec;
Y : DblDyneVec;
indx : IntDyneVec;
D, F, Prob, VarY, SDY, MeanY : double;
value, TOL, VIF, AdjR2 : double;
SSY, SSres, resvar, SSreg : double;
ColLabels : StrDyneVec;
title : string;
deplabel : string;
errcode : boolean;
begin
SetLength(X,NoCases,NoIndep+1); // augmented independent var. matrix
SetLength(XT,Noindep+1,NoCases); // transpose of independent var's
SetLength(XTX,Noindep+1,Noindep+1); // product of transpose X times X
SetLength(Y,NCases+1); // Y variable values
SetLength(XTY,Noindep+1); // X transpose time Y
SetLength(indx,Noindep+1);
SetLength(ColLabels,NCases);
// initialize
for i := 0 to Noindep do
begin
indx[i] := 0;
XTY[i] := 0.0;
Y[i] := 0.0;
tprobs[i] := 0.0;
Means[i] := 0.0;
Variances[i] := 0.0;
StdDevs[i] := 0.0;
BWeights[i] := 0.0;
BetaWeights[i] := 0.0;
errcode := false;
for j := 0 to Noindep do XTX[i,j] := 0.0;
end;
SSY := 0.0;
VarY := 0.0;
SDY := 0.0;
MeanY := 0.0;
for i := 0 to NCases-1 do
begin
ColLabels[i] := 'Case ' + IntToStr(i+1);
Y[i] := 0.0;
end;
// get independent matrix and Y vector from the grid
NCases := 0;
N := Noindep + 1;
for i := 1 to NoCases do
begin
if Not GoodRecord(i,Noindep,IndepCols) then continue;
for j := 0 to Noindep-1 do
begin
value := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[IndepCols[j],i]));
X[NCases,j] := value;
Means[j] := Means[j] + value;
Variances[j] := Variances[j] + (value * value);
end;
value := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[DepCol,i]));
Y[NCases] := value;
MeanY := MeanY + value;
SSY := SSY + (value * value);
Means[Noindep] := Means[Noindep] + value;
Variances[Noindep] := Variances[Noindep] + (value * value);
NCases := NCases + 1;
end;
deplabel := OS3MainFrm.DataGrid.Cells[DepCol,0];
RowLabels[NoIndep] := 'Intercept';
VarY := SSY - (MeanY * MeanY / NCases);
VarY := VarY / (NCases - 1);
SDY := sqrt(VarY);
title := format('Variance Y = %10.3f SSY = %10.3f SDY = %10.3f',[VarY,SSY,SDY]);
OutPutFrm.RichEdit.Lines.Add(title);
OutPutFrm.ShowModal ;
// augment the matrix
for i := 1 to NCases do X[i-1,Noindep] := 1.0;
Y[NCases] := 1.0;
// get transpose of augmented X matrix
MATTRN(XT,X,NCases,N);
{ if PrintAll then
begin
title := 'XT MATRIX';
MAT_PRINT(XT,Noindep+1,NCases,title,RowLabels,ColLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
}
// get product of the augmented X transpose times augmented X
MATAXB(XTX,XT,X,N,NCases,NCases,N,errorcode);
{ if PrintAll then
begin
title := 'XTX MATRIX';
MAT_PRINT(XTX,Noindep+1,Noindep+1,title,RowLabels,RowLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
}
//Get means, variances and standard deviations
errorcode := false;
for i := 0 to N-1 do
begin
Variances[i] := XTX[i,i] - (sqr(XTX[i,N-1]) / NCases);
Variances[i] := Variances[i] / (NCases - 1);
if (Variances[i] > 0.0) then StdDevs[i] := sqrt(Variances[i])
else errorcode := true;
Means[i] := XTX[N-1,i] / NCases;
end;
if printall then
begin
DynVectorPrint(Means,Noindep+1,'MEANS',RowLabels,NCases);
DynVectorPrint(Variances,Noindep+1,'VARIANCES',RowLabels,NCases);
DynVectorPrint(StdDevs,Noindep+1,'STD. DEV.S',RowLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
// get product of the augmented X transpose matrix times the Y vector
for i := 0 to N-1 do
begin
for j := 0 to NCases-1 do
begin
XTY[i] := XTY[i] + (XT[i,j] * Y[j]);
end;
end;
{ if PrintAll then
begin
title := 'XTY VECTOR';
DynVectorPrint(XTY,Noindep+1,title,RowLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
}
// get inverse of the augmented cross products matrix among independent variables
SVDInverse(XTX,N);
{ if PrintAll then
begin
title := 'XTX MATRIX INVERSE';
MAT_PRINT(XTX,Noindep+1,Noindep+1,title,RowLabels,RowLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
}
// multiply augmented inverse matrix times the XTY vector
// result is bweights with the intercept last
for i := 0 to N-1 do
for j := 0 to N-1 do
BWeights[i] := BWeights[i] + (XTX[i,j] * XTY[j]);
//Get Beta weightw
for i := 0 to N-2 do
BetaWeights[i] := BWeights[i] * StdDevs[i] / SDY;
// Get standard errors, squared multiple correlation, tests of significance
SSres := 0.0;
for i := 0 to NoIndep do
begin
SSres := SSres + (BWeights[i] * XTY[i]);
end;
SSres := SSY - SSres;
resvar := SSres / (NCases - N);
if resvar > 0.0 then stderrest := sqrt(resvar);
for i := 0 to N-1 do // Standard errors and t-tedt values for weights
begin
BStdErrs[i] := sqrt(resvar * XTX[i,i]);
Bttests[i] := BWeights[i] / BStdErrs[i];
tprobs[i] := probt(Bttests[i],NCases-N);
end;
SSY := VarY * (NCases-1);
SSreg := SSY - SSres;
R2 := SSreg / SSY;
F := (SSreg / (N - 1)) / (SSres / (NCases - N));
Prob := probf(F,(N-1),(NCases-N));
AdjR2 := 1.0 - (1.0 - R2) * (NCases - 1) / (NCases - N);
{
if PrintAll then
begin
OutPutFrm.RichEdit.Lines.Add('Dependent variable: ' + deplabel);
OutPutFrm.RichEdit.Lines.Add('');
title := 'B WEIGHTS';
DynVectorPrint(BWeights,Noindep+1,title,RowLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
if PrintAll then
begin
OutPutFrm.RichEdit.Lines.Add('Dependent variable: ' + deplabel);
OutPutFrm.RichEdit.Lines.Add('');
title := 'BETA WEIGHTS';
DynVectorPrint(BetaWeights,Noindep,title,RowLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
if PrintAll then
begin
title := 'B STD.ERRORS';
DynVectorPrint(BStdErrs,Noindep+1,title,RowLabels,NCases);
title := 'B t-test VALUES';
DynVectorPrint(Bttests,Noindep+1,title,RowLabels,NCases);
title := 'B t VALUE PROBABILITIES';
DynVectorPrint(tprobs,Noindep+1,title,RowLabels,NCases);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
}
if PrintAll then
begin
title := format('SSY = %10.2f, SSreg = %10.2f, SSres = %10.2f',
[SSY,SSreg,SSres]);
OutPutFrm.RichEdit.Lines.Add(title);
title := format('R2 = %6.4f, F = %8.2f, D.F. = %d %d, Prob>F = %6.4f',
[R2,F,N-1,NCases-N,Prob]);
OutPutFrm.RichEdit.Lines.Add(title);
title := format('Standard Error of Estimate = %8.2f',[stderrest]);
OutPutFrm.RichEdit.Lines.Add(title);
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
end;
RowLabels[N-1] := 'Intercept';
OutPutFrm.RichEdit.Lines.Add(' Variable Beta B Std.Err. t prob VIF TOL');
Correlations(NoIndep,IndepCols,XTX,Means,Variances, StdDevs,errcode,NCases);
SVDinverse(XTX,NoIndep);
for i := 0 to NoIndep do
begin
VIF := XTX[i,i];
if VIF > 0.0 then TOL := 1.0 / VIF else TOL := 0.0;
title := format('%10s%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f',
[RowLabels[i],BetaWeights[i],BWeights[i],BStdErrs[i],Bttests[i],
tprobs[i], VIF, TOL]);
OutPutFrm.RichEdit.Lines.Add(title);
end;
OutPutFrm.RichEdit.Lines.Add(' ');
OutPutFrm.RichEdit.Lines.Add('SOURCE DF SS MS F Prob.>F');
title := format('Regression%3d %9.3f %9.3f %6.4f',[N-1,SSreg,SSreg/(N-1),F,Prob]);
OutPutFrm.RichEdit.Lines.Add(title);
title := format('Residual %3d %9.3f %9.3f',[(NCases-1),SSres,SSres/(NCases-N)]);
OutPutFrm.RichEdit.Lines.Add(title);
title := format('Total %3d %9.3f',[NCases-1,SSY]);
OutPutFrm.RichEdit.Lines.Add(title);
OutPutFrm.RichEdit.Lines.Add('');
title := format('R2 = %6.4f, F = %8.2f, D.F. = %d %d Prob.>F = %6.4f',[R2,F,N-1,NCases-N,Prob]);
OutPutFrm.RichEdit.Lines.Add(title);
title := format('Adjusted R2 = %6.4f',[AdjR2]);
OutPutFrm.RichEdit.Lines.Add(title);
title := format('Standard Error of Estimate = %8.2f',[stderrest]);
OutPutFrm.RichEdit.Lines.Add(title);
OutPutFrm.ShowModal;
// clean up the heap
ColLabels := nil;
indx := nil;
XTY := nil;
Y := nil;
XTX := nil;
XT := nil;
X := nil;
end;
//-------------------------------------------------------------------
procedure Dynnonsymroots(var a : DblDyneMat; nv : integer;
var nf : integer; c : real;
var v : DblDyneMat; var e : DblDyneVec;
var px : DblDyneVec;
var t : double;
var ev : double);
{ roots and vectors of a non symetric matrix. a is square matrix entered
and is destroyed in process. nv is number of variables (rows and columns )
of a. nf is the number of factorsto be extracted - is output as the number
which exceeded c, the minimum eigenvalue to be extracted. v is the output
matrix of column vectors of loadings. e is the output vector of roots. px
is the percentages of trace for factors. t is the trace of the matrix and
ev is the percent of trace extracted }
label 40;
var
y, z : DblDyneVec;
ek, e2, d : real;
i, j, k, m : integer;
begin
SetLength(y,nv);
SetLength(z,nv);
t := 0.0;
for i := 0 to nv-1 do t := t + a[i,i];
for k := 0 to nf-1 do
begin
for i := 0 to nv-1 do
begin
px[i] := 1.0;
y[i] := 1.0;
end;
e[k] := 1.0;
ek := 1.0;
for m := 1 to 25 do
begin
for i := 0 to nv-1 do
begin
v[i,k] := px[i] / e[k];
z[i] := y[i] / ek;
end;
for i := 0 to nv - 1 do
begin
px[i] := 0.0;
for j := 0 to nv-1 do px[i] := px[i] + a[i,j] * v[j,k];
y[i] := 0.0;
for j := 0 to nv-1 do y[i] := y[i] + a[j,i] * z[j];
end;
e2 := 0.0;
for j := 0 to nv-1 do e2 := e2 + px[j] * v[j,k];
e[k] := sqrt(abs(e2));
ek := 0.0;
for j := 0 to nv-1 do ek := ek + y[j] * z[j];
ek := sqrt(abs(ek));
end;
if e2 >= sqr(c) then
begin
d := 0.0;
for j := 0 to nv - 1 do d := d + v[j,k] * z[j];
d := e[k] / d;
for i := 0 to nv - 1 do
for j := 0 to nv - 1 do
a[i,j] := a[i,j] - v[i,k] * z[j] * d;
end
else begin
nf := k - 1;
goto 40;
end;
end;
40 : for i := 0 to nf-1 do px[i] := e[i] / t * 100.0;
ev := 0.0;
for i := 0 to nf-1 do ev := ev + px[i];
z := nil;
y := nil;
end; { of procedure nonsymroots }
//-----------------------------------------------------------------------------
function DynCorrelations(novars : integer;
VAR ColSelected : IntDyneVec;
VAR DataGrid : DblDyneMat;
VAR rmatrix : DblDyneMat;
VAR means : DblDyneVec;
VAR vars : DblDyneVec;
VAR stddevs : DblDyneVec;
NCases : integer;
ReturnType : integer) : integer;
var
i, j, k, row, col, errorcode : integer;
X, Y : double;
begin
errorcode := 0;
for i := 0 to novars - 1 do
begin
means[i] := 0.0;
vars[i] := 0.0;
stdDevs[i] := 0.0;
for j := 0 to novars - 1 do
begin
rmatrix[i,j] := 0.0;
end;
end;
{ get cross products }
for i := 0 to NCases - 1 do
begin
if IsFiltered(i) then continue;
for j := 0 to novars - 1 do
begin
row := ColSelected[j];
X := DataGrid[i,row];
means[j] := means[j] + X;
vars[j] := vars[j] + (X * X);
for k := 0 to novars - 1 do
begin
col := ColSelected[k];
Y := DataGrid[i,col];
rmatrix[j,k] := rmatrix[j,k] + (X * Y);
end;
end;
end;
for j := 0 to novars - 1 do
begin
vars[j] := vars[j] - (means[j] * means[j] / NCases);
vars[j] := vars[j] / (NCases-1);
if (vars[j] > 0.0) then stddevs[j] := sqrt(vars[j])
else stddevs[j] := 0.0;
end;
if ReturnType = 1 then {return cross-products, variances, std.devs, means }
begin
for i := 0 to novars - 1 do
begin
means[i] := means[i] / NCases;
end;
DynCorrelations := errorcode;
exit;
end;
for i := 0 to novars - 1 do {get variance-covariance matrix }
begin
for j := 0 to novars - 1 do
begin
rmatrix[i,j] := rmatrix[i,j] - (means[i] * means[j] / NCases);
rmatrix[i,j] := rmatrix[i,j] / (NCases - 1);
end;
end;
if ReturnType = 2 then
begin
for i := 0 to novars - 1 do
begin
means[i] := means[i] / NCases;
end;
DynCorrelations := errorcode;
exit;
end;
for i := 0 to novars - 1 do { get product-moment correlations }
begin
for j := 0 to novars - 1 do
begin
if ((stddevs[i] > 0.0) and (stddevs[j] > 0.0)) then
rmatrix[i,j] := rmatrix[i,j] / (stddevs[i] * stddevs[j])
else
begin
rmatrix[i,j] := 9.999;
errorcode := 1;
end;
end;
end;
for i := 0 to novars - 1 do
begin
means[i] := means[i] / NCases;
end;
DynCorrelations := errorcode;
end;
//---------------------------------------------------------------------------
procedure Predict(VAR ColNoSelected : IntDyneVec;
NoVars : integer;
VAR IndepInverse : DblDyneMat;
VAR Means : DblDyneVec;
VAR StdDevs : DblDyneVec;
VAR BetaWeights : DblDyneVec;
StdErrEst : double;
VAR IndepIndex : IntDyneVec;
NoIndepVars : integer);
var
col, i, j, k, index, IndexX, IndexY : integer;
predicted, zpredicted, z1, z2, resid, Term1, Term2 : double;
StdErrPredict, t95, Hi95, Low95 : double;
begin
{ use the next available grid column to store the z predicted score }
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'Pred.z';
OS3MainFrm.DataGrid.Cells[col,0] := 'Pred.z';
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'z Resid.';
OS3MainFrm.DataGrid.Cells[col,0] := 'z Resid.';
for i := 1 to NoCases do
begin
zpredicted := 0.0;
for j := 1 to NoIndepVars do
begin
Index := IndepIndex[j-1];
k := ColNoSelected[Index-1];
z1 := (StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[k,i])) -
Means[Index-1]) / StdDevs[index-1];
zpredicted := zpredicted + (z1 * BetaWeights[j-1]);
OS3MainFrm.DataGrid.Cells[col-1,i] := format('%8.4f',[zpredicted]);
end;
Index := ColNoSelected[NoVars-1];
z2 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[Index,i]));
z2 := (z2 - Means[NoVars-1]) / StdDevs[NoVars-1];
OS3MainFrm.DataGrid.Cells[col,i] := format('%8.4f',[(z2 - zpredicted)]);
end;
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'Pred.Raw';
OS3MainFrm.DataGrid.Cells[col,0] := 'Pred.Raw';
{ calculate raw predicted scores and store in grid at col }
for i := 1 to NoCases do
begin
predicted := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[col-2,i])) *
StdDevs[NoVars-1] + Means[NoVars-1];
OS3MainFrm.DataGrid.Cells[col,i] := format('%8.3f',[predicted]);
end;
{ Calculate residuals of predicted raw scores }
col := NoVariables +1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'Raw Resid.';
OS3MainFrm.DataGrid.Cells[col,0] := 'Raw Resid.';
for i := 1 to NoCases do
begin
Index := ColNoSelected[NoVars-1];
resid := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[col-1,i])) -
StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[Index,i]));
OS3MainFrm.DataGrid.Cells[col,i] := Format('%8.3f',[resid]);
end;
{ Calculate Confidence Interval for raw predicted score }
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'StdErrPred';
OS3MainFrm.DataGrid.Cells[col,0] := 'StdErrPred';
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'Low 95%';
OS3MainFrm.DataGrid.Cells[col,0] := 'Low 95%';
col := NoVariables + 1;
DictionaryFrm.NewVar(col);
DictionaryFrm.DictGrid.Cells[1,col] := 'Top 95%';
OS3MainFrm.DataGrid.Cells[col,0] := 'Top 95%';
for i := 1 to NoCases do
begin
{ get term1 of the std. err. prediction }
Term1 := 0.0;
for j := 1 to NoIndepVars do
begin
Index := IndepIndex[j-1];
col := ColNoSelected[Index-1];
z1 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[col,i]));
z1 := (z1 - Means[Index-1]) / StdDevs[Index-1];
z1 := (z1 * z1) * IndepInverse[j-1,j-1];
Term1 := Term1 + z1;
end;
{ get term2 of the std err. of prediction }
term2 := 0.0;
for j := 1 to NoIndepVars - 1 do
begin
for k := j + 1 to NoIndepVars do
begin
IndexX := IndepIndex[j-1];
col := ColNoSelected[IndexX-1];
z1 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[col,i]));
IndexY := IndepIndex[k-1];
col := ColNoSelected[IndexY-1];
z2 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[col,i]));
z1 := (z1 - Means[IndexX-1]) / StdDevs[IndexX-1];
z2 := (z2 - Means[IndexY-1]) / StdDevs[IndexY-1];
Term2 := Term2 + IndepInverse[j-1,k-1] * z1 * z2;
end;
end;
term2 := 2.0 * Term2;
StdErrPredict := sqrt(NoCases + 1 + Term1 + Term2);
StdErrPredict := (StdErrEst / sqrt(NoCases)) * StdErrPredict;
t95 := Inverset(0.975,NoCases-NoIndepVars-1);
low95 := StrToFloat(Trim(OS3MainFrm.DataGrid.Cells[NoVars+4,i]));
hi95 := low95;
low95 := low95 - (t95 * StdErrPredict);
hi95 := hi95 + (t95 * StdErrPredict);
OS3MainFrm.DataGrid.Cells[NoVariables,i] := Format('%8.3f',[hi95]);
OS3MainFrm.DataGrid.Cells[NoVariables-1,i] := Format('%8.3f',[low95]);
OS3MainFrm.DataGrid.Cells[NoVariables-2,i] := Format('%8.3f',[StdErrPredict]);
end; { next case }
end;
//---------------------------------------------------------------------------
procedure MReg2(NCases : integer;
NoVars : integer;
VAR NoIndepVars : integer;
VAR IndepIndex : IntDyneVec;
VAR corrs : DblDyneMat;
VAR IndepCorrs : DblDyneMat;
VAR RowLabels : StrDyneVec;
VAR R2 : double;
VAR BetaWeights : DblDyneVec;
VAR Means : DblDyneVec;
VAR Variances : DblDyneVec;
VAR errorcode : integer;
VAR StdErrEst : double;
VAR constant : double;
probout : double;
Printit : boolean;
TestOut : boolean;
PrintInv : boolean);
{
The following routine obtains multiple regression results for a
correlation matrix consisting of 1 to NoVars. The last variable
represents the dependent variable. The number of independent
variables is passed as NoIndepVars. The inverse matrix of independent
variables may be obtained by the calling program using the variable
IndepCorrs. The user may request printing of the inverse using the
boolean variable Printit.
}
var
i, j, k, l : integer;
IndexX, IndexY : integer;
IndRowLabels : StrDyneVec;
IndColLabels : StrDyneVec;
XYCorrs : DblDyneVec;
df1, df2, df3 : double;
SSt, SSres, SSreg : double;
VarEst, F : double;
FprobF : double;
outline : string;
valstring : string;
title : string;
deplabel : string;
sum, B, Beta : double;
SSx, StdErrB : double;
AdjR2 : double;
VIF, TOL : double;
outcount : integer;
varsout : IntDyneVec;
begin
SetLength(IndRowLabels,NoVars);
SetLength(IndColLabels,NoVars);
SetLength(XYCorrs,NoVars);
SetLength(varsout,NoVars);
errorcode := 0;
outcount := 0;
VIF := 0.0;
deplabel := RowLabels[NoVars-1];
for i := 0 to NoIndepVars-1 do
begin
IndexX := IndepIndex[i];
for j := 0 to NoIndepVars-1 do
begin
IndexY := IndepIndex[j];
IndepCorrs[i,j] := corrs[IndexX-1,IndexY-1];
end;
end;
for i := 0 to NoIndepVars-1 do
begin
IndRowLabels[i] := RowLabels[IndepIndex[i]-1];
IndColLabels[i] := RowLabels[IndepIndex[i]-1];
XYCorrs[i] := corrs[IndepIndex[i]-1,NoVars-1];
end;
SVDinverse(IndepCorrs,NoIndepVars);
title := 'Inverse of independent variables matrix';
if (PrintInv = true) then
MAT_PRINT(IndepCorrs,NoIndepVars,NoIndepVars,title,IndRowLabels,IndColLabels,NCases);
{ Get product of inverse matrix times vector of correlations
between independent and dependent variables }
R2 := 0.0;
for i := 0 to NoIndepVars-1 do
begin
BetaWeights[i] := 0.0;
for j := 0 to NoIndepVars-1 do
begin
BetaWeights[i] := BetaWeights[i] + IndepCorrs[i,j] * XYCorrs[j];
end;
R2 := R2 + BetaWeights[i] * XYCorrs[i];
end;
df1 := NoIndepVars;
df2 := NCases - NoIndepVars - 1;
df3 := NCases - 1;
SSt := (NCases-1) * Variances[NoVars-1];
SSres := SSt * (1.0 - R2);
SSreg := SSt - SSres;
VarEst := SSres / df2;
if (VarEst > 0.0) then StdErrEst := sqrt(VarEst)
else
begin
ShowMessage('ERROR! Error in computing variance estimate.');
StdErrEst := 0.0;
end;
if (R2 < 1.0) and (df2 > 0.0) and (df1 > 0.0) then F := (R2 / df1) / ((1.0-R2)/ df2)
else F := 0.0;
FProbF := probf(F,df1,df2);
OutPutFrm.RichEdit.Lines.Add('SOURCE DF SS MS F Prob.>F');
outline := format('Regression %3.0f %9.3f %9.3f %9.3f %9.3f',[df1,SSreg,SSreg/df1,F,FprobF]);
OutPutFrm.RichEdit.Lines.Add(outline);
outline := format('Residual %3.0f %9.3f %9.3f',[df2,SSres,SSres/df2]);
OutPutFrm.RichEdit.Lines.Add(outline);
outline := format('Total %3.0f %9.3f',[df3,SSt]);
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('');
AdjR2 := 1.0 - (1.0 - R2) * (NCases - 1) / df2;
if (PrintIt = true) then
begin
// OutPutFrm.RichEdit.ParaGraph.Alignment := taLeftJustify;
OutPutFrm.RichEdit.Lines.Add('Dependent Variable: ' + deplabel);
OutPutFrm.RichEdit.Lines.Add('');
outline := format('%8s%10s%10s%12s%5s%5s',['R','R2','F','Prob.>F','DF1','DF2']);
OutPutFrm.RichEdit.Lines.Add(outline);
outline := format('%8.3f%10.3f%10.3f%10.3f%5.0f%5.0f',
[sqrt(R2),R2,F,FProbF,df1,df2]);
OutPutFrm.RichEdit.Lines.Add(outline);
outline := format('Adjusted R Squared = %5.3f',[AdjR2]);
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('');
outline := format('Std. Error of Estimate = %10.3f',[StdErrEst]);
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add('Variable Beta B Std.Error t Prob.>t VIF TOL');
end;
df1 := 1.0;
df2 := NCases - NoIndepVars - 1;
sum := 0.0;
for i := 0 to NoIndepVars-1 do
begin
beta := BetaWeights[i];
B := beta * sqrt(Variances[NoVars-1]) / sqrt(Variances[IndepIndex[i]-1]);
sum := sum + B * Means[IndepIndex[i]-1];
SSx := (NCases-1) * Variances[IndepIndex[i]-1];
if (IndepCorrs[i,i] > 0.0) and (VarEst > 0.0) then
begin
StdErrB := sqrt(VarEst / (SSx * (1.0 / IndepCorrs[i,i])));
F := B / StdErrB;
FProbF := probf(F*F,df1,df2);
VIF := IndepCorrs[i,i];
TOL := 1.0 / VIF;
end
else
begin
ShowMessage('ERROR! Error in estimating std.err. of a B');
StdErrB := 0.0;
F := 0.0;
FProbF := 0.0;
end;
if (PrintIt = true) then
begin
valstring := format('%10s',[IndRowLabels[i]]);
outline := format('%10s%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f',
[valstring, beta ,B, StdErrB, F, FProbF, VIF, TOL]);
if FprobF > ProbOut then outline := outline + ' Exceeds limit - to be removed.';
OutPutFrm.RichEdit.Lines.Add(outline);
end;
if FprobF > ProbOut then
begin
outcount := outcount + 1;
varsout[outcount-1] := IndepIndex[i];
end;
end;
if (PrintIt = true) then OutPutFrm.RichEdit.Lines.Add('');
{ Get constant }
constant := Means[NoVars-1] - sum;
if (PrintIt = true) then
begin
outline := format('Constant = %10.3f',[constant]);
OutPutFrm.RichEdit.Lines.Add(outline);
end;
{ Now remove any variables that exceed tolerance }
if (outcount > 0) and (TestOut = true) then
begin
for i := 0 to outcount-1 do
begin
k := varsout[i]; { variable to eliminate }
for j := 0 to NoIndepVars-1 do
begin
if IndepIndex[j] = k then {eliminate this one }
begin
for l := j to NoIndepVars-2 do
IndepIndex[l] := IndepIndex[l+1];
end;
end;
end;
NoIndepVars := NoIndepVars - outcount;
errorcode := outcount;
end;
varsout := nil;
XYCorrs := nil;
IndColLabels := nil;
IndRowLabels := nil;
end;
//---------------------------------------------------------------------------
procedure MATSUB(VAR a, b, c : DblDyneMat;
brows, bcols, crows, ccols : integer; VAR errorcode : boolean);
// Subtracts matrix c from b and returns the results in matrix a
var i, j : integer;
begin
errorcode := FALSE;
if ((brows <> crows) or (bcols <> ccols)) then errorcode := TRUE
else
begin
for i := 0 to brows-1 do
for j := 0 to bcols-1 do
a[i,j] := b[i,j] - c[i,j];
end;
end; { of matsub }
//---------------------------------------------------------------------------
procedure IntArrayPrint(mat : IntDyneMat;
rows, cols : integer;
ytitle : string;
RowLabels, ColLabels : StrDyneVec;
Title : string);
var
i, j, first, last, nflds : integer;
done : boolean;
outline: string;
valstring: string;
begin
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add(Title);
OutPutFrm.RichEdit.Lines.Add('');
nflds := 4;
done := FALSE;
first := 1;
while not done do
begin
OutPutFrm.RichEdit.Lines.Add('');
outline := ' ' + ytitle;
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('Variables');
outline := ' ';
last := first + nflds;
if last >= cols then
begin
done := TRUE;
last := cols
end;
for i := first to last do
begin
outline := outline + format('%13s',[ColLabels[i-1]]);
end;
OutPutFrm.RichEdit.Lines.Add(outline);
for i := 1 to rows do
begin
outline := format('%10s',[RowLabels[i-1]]);
for j := first to last do
begin
valstring := format('%12d ',[mat[i-1,j-1]]);
outline := outline + valstring;
end;
OutPutFrm.RichEdit.Lines.Add(outline);
end;
OutPutFrm.RichEdit.Lines.Add('');
first := last + 1
end;
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add('');
end;
//---------------------------------------------------------------------------
procedure eigens(VAR a: DblDyneMat; Var d : DblDyneVec; n : integer);
var e : DblDyneVec;
i : integer;
begin
SetLength(e,n);
for i := 1 to n do
begin
d[i-1] := 0.0;
e[i-1] := 0.0;
end;
tred2(a, n, d ,e ); { Upon return, d contains diagonal values, e contains
off diagonal values, and a contains the orthogonal
matrix from the tridiagonalization of the a matrix }
tqli(d, e, n, a); { Upon return, d contains eigenvalues, a the column
eigenvectors of the tridiagonal matrix (in d and e). }
e := nil;
end; { Procedure eigens }
//-------------------------------------------------------------------
PROCEDURE tred2(VAR a: DblDyneMat; n: integer; VAR d,e: DblDyneVec);
(* Programs using routine TRED2 must define the types
TYPE
glnp = ARRAY [1..np] OF real;
glnpnp = ARRAY [1..np,1..np] OF real;
where 'np by np' is the physical dimension of the matrix to be analyzed. *)
VAR
l,k,j,i: integer;
scale,hh,h,g,f: double;
BEGIN
IF (n > 1) THEN BEGIN
FOR i := n DOWNTO 2 DO BEGIN
l := i-1;
h := 0.0;
scale := 0.0;
IF (l > 1) THEN BEGIN
FOR k := 1 to l DO BEGIN
scale := scale+abs(a[i-1,k-1])
END;
IF (scale = 0.0) THEN BEGIN
e[i-1] := a[i-1,l-1]
END ELSE BEGIN
FOR k := 1 to l DO BEGIN
a[i-1,k-1] := a[i-1,k-1]/scale;
h := h+sqr(a[i-1,k-1])
END;
f := a[i-1,l-1];
g := -sign(sqrt(h),f);
e[i-1] := scale*g;
h := h-f*g;
a[i-1,l-1] := f-g;
f := 0.0;
FOR j := 1 to l DO BEGIN
(* Next statement can be omitted if eigenvectors not wanted *)
a[j-1,i-1] := a[i-1,j-1]/h;
g := 0.0;
FOR k := 1 to j DO BEGIN
g := g+a[j-1,k-1]*a[i-1,k-1]
END;
IF (l > j) THEN FOR k := j+1 to l DO g := g+a[k-1,j-1]*a[i-1,k-1];
e[j-1] := g/h;
f := f+e[j-1]*a[i-1,j-1]
END;
hh := f/(h+h);
FOR j := 1 to l DO BEGIN
f := a[i-1,j-1];
g := e[j-1]-hh*f;
e[j-1] := g;
FOR k := 1 to j DO a[j-1,k-1] := a[j-1,k-1]-f*e[k-1]-g*a[i-1,k-1]
END
END
END ELSE BEGIN
e[i-1] := a[i-1,l-1]
END;
d[i-1] := h
END
END;
(* Next statement can be omitted if eigenvectors not wanted *)
d[0] := 0.0;
e[0] := 0.0;
FOR i := 1 to n DO BEGIN
(* Contents of this loop can be omitted if eigenvectors not wanted,
except for statement d[i] := a[i,i]; *)
l := i-1;
IF (d[i-1] <> 0.0) THEN BEGIN
FOR j := 1 to l DO BEGIN
g := 0.0;
FOR k := 1 to l DO BEGIN
g := g+a[i-1,k-1]*a[k-1,j-1]
END;
FOR k := 1 to l DO BEGIN
a[k-1,j-1] := a[k-1,j-1]-g*a[k-1,i-1]
END
END
END;
d[i-1] := a[i-1,i-1];
a[i-1,i-1] := 1.0;
IF (l >= 1) THEN BEGIN
FOR j := 1 to l DO BEGIN
a[i-1,j-1] := 0.0;
a[j-1,i-1] := 0.0
END
END
END
END;
//-------------------------------------------------------------------
PROCEDURE tqli(VAR d,e: DblDyneVec; n: integer; VAR z: DblDyneMat);
LABEL 1,2;
VAR
m,l,iter,i,k: integer;
s,r,p,g,f,dd,c,b: double;
BEGIN
IF (n > 1) THEN BEGIN
FOR i := 2 to n DO BEGIN
e[i-2] := e[i-1]
END;
e[n-1] := 0.0;
FOR l := 1 to n DO BEGIN
iter := 0;
1: FOR m := l to n-1 DO BEGIN
dd := abs(d[m-1])+abs(d[m]);
IF (abs(e[m-1])+ dd = dd) THEN GOTO 2
END;
m := n;
2: IF (m <> l) THEN BEGIN
IF (iter = 30) THEN BEGIN
ShowMessage('Too many iterations in routine tqli');
exit;
END;
iter := iter+1;
g := (d[l]-d[l-1])/(2.0*e[l-1]);
r := sqrt(sqr(g)+1.0);
g := d[m-1] - d[l-1] + e[l-1] / (g+sign(r,g));
s := 1.0;
c := 1.0;
p := 0.0;
FOR i := m-1 DOWNTO l DO BEGIN
f := s * e[i-1];
b := c * e[i-1];
IF (abs(f) >= abs(g)) THEN BEGIN
c := g / f;
r := sqrt(sqr(c) + 1.0);
e[i] := f * r;
s := 1.0 / r;
c := c * s
END ELSE BEGIN
s := f / g;
r := sqrt(sqr(s) + 1.0);
e[i] := g * r;
c := 1.0 / r;
s := s * c
END;
g := d[i] - p;
r := (d[i-1] - g) * s + 2.0 * c * b;
p := s * r;
d[i] := g + p;
g := c * r - b;
(* Next loop can be omitted if eigenvectors not wanted *)
FOR k := 1 to n DO BEGIN
f := z[k-1,i];
z[k-1,i] := s * z[k-1,i-1] + c * f;
z[k-1,i-1] := c * z[k-1,i-1] - s * f
END
END;
d[l-1] := d[l-1] - p;
e[l-1] := g;
e[m-1] := 0.0;
GOTO 1
END
END
END
END;
//-------------------------------------------------------------------
function SEVS(nv,nf : integer;
c : double;
var r : DblDyneMat;
VAR v : DblDyneMat;
VAR e : DblDyneVec;
var p : DblDyneVec;
VAR nd : integer) : integer ;
{ extracts roots and denormal vectors from a symetric matrix. Veldman, 1967,
page 209 }
label 1,2;
var t, ee, ev : double;
i, j, k, m : integer;
begin
t := 0.0;
for i := 1 to nv do t := t + r[i-1,i-1];
for k := 1 to nf do { compute roots in e[k] and vector in v^[.k] }
begin
for i := 1 to nv do p[i-1] := 1.0;
begin
e[k-1] := 1.0;
for m := 1 to 25 do
begin
for i := 1 to nv do v[i-1,k-1] := p[i-1] / e[k-1];
for i := 1 to nv do p[i-1] := SCPF(r,v,-i,k,nv,nd);
ee := 0.0;
for j := 1 to nv do ee := ee + p[j-1] * v[j-1,k-1];
e[k-1] := sqrt(abs(ee));
end;
end;
if ee < (c * c) then goto 1;
for i := 1 to nv do
for j := 1 to nv do
r[i-1,j-1] := r[i-1,j-1] - (v[i-1,k-1] * v[j-1,k-1]);
end;
goto 2;
1 : nf := k - 1;
2 : for i := 1 to nf do p[i-1] := e[i-1] / t * 100.0;
ev := 0.0;
for i := 1 to nf do ev := ev + p[i-1];
{
stopit;
writeln(lst);
writeln(lst,'Root % Extracted');
for i := 1 to nf do writeln(lst,i:3,' ',p^[i]:6:3);
writeln(lst,' Trace = ',t:6:3,' % Extracted = ',ev:6:3);
writeln(lst);
}
result := nf;
end; { of SEVS procedure }
//-------------------------------------------------------------------
function SCPF(VAR x,y : DblDyneMat; kx,ky,n,nd : integer) : double;
{ sum of cross products of two vectors. Veldman, 1967, pp 128-129 }
var j,k,i : integer;
scp : double;
begin
scp := 0.0;
scpf := 0.0;
j := abs(kx);
k := abs(ky);
if ((kx = 0) and (ky = 0)) then exit;
if ((kx < 0) and (ky < 0)) then
begin
for i := 1 to n do scp := scp + x[j-1,i-1] * y[k-1,i-1];
end;
if ((kx < 0) and (ky > 0)) then
begin
for i := 1 to n do scp := scp + x[j-1,i-1] * y[i-1,k-1];
end;
if ((kx > 0) and (ky < 0)) then
begin
for i := 1 to n do scp := scp + x[i-1,j-1] * y[k-1,i-1];
end;
if ((kx > 0) and (ky > 0)) then
begin
for i := 1 to n do scp := scp + x[i-1,j-1] * y[i-1,k-1];
end;
scpf := scp;
end; { of SCPF }
//-------------------------------------------------------------------
procedure MAT_PRINT(VAR xmat : DblDyneMat ;
rows, cols : integer;
VAR title : string;
VAR RowLabels: StrDyneVec;
VAR ColLabels: StrDyneVec;
NCases: integer);
var
i, j, first, last, nflds : integer;
done : boolean;
outline: string;
valstring: string;
begin
OutPutFrm.RichEdit.Lines.Add('');
outline := format(' with %4d cases.',[NCases]);
outline := title + outline;
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('');
nflds := 4;
done := FALSE;
first := 1;
while not done do
begin
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add('Variables');
outline := ' ';
last := first + nflds;
if last >= cols then
begin
done := TRUE;
last := cols
end;
for i := first to last do
begin
outline := outline + format('%13s',[ColLabels[i-1]]);
end;
OutPutFrm.RichEdit.Lines.Add(outline);
for i := 1 to rows do
begin
outline := format('%10s',[RowLabels[i-1]]);
for j := first to last do
begin
valstring := format('%12.3f ',[xmat[i-1,j-1]]);
outline := outline + valstring;
end;
OutPutFrm.RichEdit.Lines.Add(outline);
end;
OutPutFrm.RichEdit.Lines.Add('');
first := last + 1
end;
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add('');
end;
//--------------------------------------------------------------------
procedure DynVectorPrint(VAR avector: DblDyneVec;
NoVars : integer;
title : string;
VAR Labels : StrDyneVec;
NCases : integer);
var
i, j, first, last, nflds : integer;
done : boolean;
outline: string;
valstring: string;
begin
OutPutFrm.RichEdit.Lines.Add('');
outline := format(' with %4d valid cases.',[NCases]);
outline := title + outline;
OutPutFrm.RichEdit.Lines.Add(outline);
nflds := 4;
done := FALSE;
first := 0;
while not done do
begin
OutPutFrm.RichEdit.Lines.Add('');
outline := 'Variables';
last := first + nflds;
if last >= NoVars -1 then
begin
done := TRUE;
last := NoVars-1;
end;
for i := first to last do
begin
outline := outline + format('%13s',[Labels[i]]);
end;
OutPutFrm.RichEdit.Lines.Add(outline);
outline := ' ';
for j := first to last do
begin
valstring := format('%12.3f ',[avector[j]]);
outline := outline + valstring;
end;
OutPutFrm.RichEdit.Lines.Add(outline);
first := last + 1
end;
OutPutFrm.RichEdit.Lines.Add('');
end;
//--------------------------------------------------------------------------
procedure scatplot(var x : DblDyneVec;
var y : DblDyneVec;
nocases : integer;
titlestr : string;
x_axis, y_axis : string;
x_min, x_max, y_min, y_max : double;
VAR VarLabels : StrDyneVec);
var
i, j, k, l, row, xslot : integer;
xdelta, ypred, xtemp, maxx, maxy, minx, miny, stepx, stepy : double;
incrementx, incrementy, rangex, rangey, swap : double;
plotstring : array[0..51,0..61] of char;
ymed, xmed : double;
height : integer;
overlap : boolean;
valuestring : string[2];
howlong : integer;
outline : string;
Labels : StrDyneVec;
begin
SetLength(Labels,NoVariables);
for i := 1 to nocases do Labels[i-1] := VarLabels[i-1];
height := 40;
rangex := x_max - x_min ;
incrementx := rangex / 15.0;
xdelta := rangex / 60;
xmed := rangex / 2;
rangey := y_max - y_min;
incrementy := rangey / height;
ymed := rangey / 2;
{ sort in descending order }
for i := 1 to (nocases - 1) do
begin
for j := (i + 1) to nocases do
begin
if y[i-1] < y[j-1] then
begin
swap := y[i-1];
y[i-1] := y[j-1];
y[j-1] := swap;
swap := x[i-1];
x[i-1] := x[j-1];
x[j-1] := swap;
outline := Labels[i-1];
Labels[i-1] := Labels[j-1];
Labels[j-1] := outline;
end;
end;
end;
outline := ' SCATTERPLOT - ' + titlestr;
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add(y_axis);
maxy := y_max;
for i := 1 to 60 do
for j := 1 to height+1 do plotstring[j,i] := ' ';
{ Set up the plot strings with the data }
row := 0;
while maxy > y_min do
begin
row := row + 1;
plotstring[row,30] := '|';
if (row = (height / 2)) then
begin
for i := 1 to 60 do plotstring[row,i] := '-';
end;
for i := 1 to nocases do
begin
if ((maxy >= y[i-1]) and (y[i-1] > (maxy - incrementy))) then
begin
xslot := round(((x[i-1] - x_min) / rangex) * 60);
if xslot < 1 then xslot := 1;
if xslot > 60 then xslot := 60;
overlap := false;
str(i:2,valuestring);
howlong := 1;
if (valuestring[1] <> ' ') then howlong := 2;
for l := xslot to (xslot + howlong - 1) do
if (plotstring[row,l] = '*') then overlap := true;
if (overlap) then plotstring[row,xslot] := '*'
else
begin
if (howlong < 2) then
plotstring[row,xslot] := valuestring[2]
else for l := 1 to 2 do
plotstring[row,xslot + l - 1] := valuestring[l];
end;
end;
end;
maxy := maxy - incrementy;
end;
{ print the plot }
for i := 1 to row do
begin
outline := ' |';
for j := 1 to 60 do outline := outline + format('%1s',[plotstring[i,j]]);
outline := outline + format('|-%6.2f-%6.2f',
[(y_max - i * incrementy),(y_max - i * incrementy + incrementy)]);
OutPutFrm.RichEdit.Lines.Add(outline);
end;
outline := '';
for i := 1 to 63 do outline := outline + '-';
OutPutFrm.RichEdit.Lines.Add(outline);
outline := '';
for i := 1 to 16 do outline := outline + ' | ';
outline := outline + x_axis;
OutPutFrm.RichEdit.Lines.Add(outline);
outline := '';
for i := 1 to 16 do outline := outline + format('%4.1f',[(x_min + i * incrementx - incrementx)]);
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add('Labels:');
for i := 1 to nocases do
begin
outline := format('%2d = %s',[i,Labels[i-1]]);
OutPutFrm.RichEdit.Lines.Add(outline);
end;
OutPutFrm.ShowModal;
OutPutFrm.RichEdit.Clear;
Labels := nil;
end; { of scatplot procedure }
//-------------------------------------------------------------------
procedure DynIntMatPrint(mat : IntDyneMat;
rows, cols : integer;
ytitle : string;
RowLabels, ColLabels : StrDyneVec;
Title : string);
var
i, j, first, last, nflds : integer;
done : boolean;
outline: string;
valstring: string;
begin
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add(Title);
OutPutFrm.RichEdit.Lines.Add('');
nflds := 4;
done := FALSE;
first := 0;
while not done do
begin
OutPutFrm.RichEdit.Lines.Add('');
outline := ' ' + ytitle;
OutPutFrm.RichEdit.Lines.Add(outline);
OutPutFrm.RichEdit.Lines.Add('Variables');
outline := ' ';
last := first + nflds;
if last >= cols-1 then
begin
done := TRUE;
last := cols-1
end;
for i := first to last do
begin
outline := outline + format('%13s',[ColLabels[i]]);
end;
OutPutFrm.RichEdit.Lines.Add(outline);
for i := 0 to rows-1 do
begin
outline := format('%10s',[RowLabels[i]]);
for j := first to last do
begin
valstring := format('%12d ',[mat[i,j]]);
outline := outline + valstring;
end;
OutPutFrm.RichEdit.Lines.Add(outline);
end;
OutPutFrm.RichEdit.Lines.Add('');
first := last + 1
end;
OutPutFrm.RichEdit.Lines.Add('');
OutPutFrm.RichEdit.Lines.Add('');
end;
procedure SymMatRoots(A: DblDyneMat; M: integer; var E: DblDyneVec;
var V: DblDyneMat);
Label one, three, nine, fifteen;
var
L, IT, j, k : integer;
Test, sum1, sum2 : double;
X, Y, Z : DblDyneVec;
begin
// Adapted from: "Multivariate Data Analysis" by William W. Cooley and Paul
// R. Lohnes, 1971, page 121
SetLength(X, M);
SetLength(Y, M);
SetLength(Z, M);
sum2 := 0.0;
L := 0;
Test := 0.00000001;
one:
IT := 0;
for j := 0 to M-1 do Y[j] := 1.0;
three:
IT := IT + 1;
for j := 0 to M-1 do
begin
X[j] := 0.0;
for k := 0 to M-1 do X[j] := X[j] + (A[j,k] * Y[k]);
end;
E[L] := X[0];
Sum1 := 0.0;
for j := 0 to M-1 do
begin
V[j,L] := X[j] / X[0];
Sum1 := Sum1 + abs(Y[j] - V[j,L]);
Y[j] := V[j,L];
end;
if (IT - 10) <> 0 then goto nine;
if (Sum2 - Sum1) > 0 then goto nine
else
begin
showmessage('Root not converging. Exiting.');
exit;
end;
nine:
Sum2 := Sum1;
if (Sum1 - Test) > 0 then goto three;
Sum1 := 0.0;
for j := 0 to M-1 do Sum1 := Sum1 + (V[j,L] * V[j,L]);
Sum1 := sqrt(Sum1);
for j := 0 to M-1 do V[j,L] := V[j,L] / Sum1;
for j := 0 to M-1 do
for k := 0 to M-1 do
A[j,k] := A[j,k] - (V[j,L] * V[k,L] * E[L]);
if ((M-1)-L) <= 0 then goto fifteen;
L := L + 1;
goto one;
fifteen:
Z := nil;
Y := nil;
X := nil;
end;
procedure matinv(a, vtimesw, v, w: DblDyneMat; n: integer);
LABEL 1,2,3;
VAR
ainverse : array of array of double;
m,mp,np,nm,l,k,j,its,i: integer;
z,y,x,scale,s,h,g,f,c,anorm: double;
rv1: array of double;
begin
setlength(rv1,n);
setlength(ainverse,n,n);
m := n;
mp := n;
np := n;
g := 0.0;
scale := 0.0;
anorm := 0.0;
FOR i := 0 to n-1 DO BEGIN
l := i+1;
rv1[i] := scale*g;
g := 0.0;
s := 0.0;
scale := 0.0;
IF (i <= m-1) THEN BEGIN
FOR k := i to m-1 DO BEGIN
scale := scale+abs(a[k,i])
END;
IF (scale <> 0.0) THEN BEGIN
FOR k := i to m-1 DO BEGIN
a[k,i] := a[k,i]/scale;
s := s+a[k,i]*a[k,i]
END;
f := a[i,i];
g := -sign(sqrt(s),f);
h := f*g-s;
a[i,i] := f-g;
IF (i <> n-1) THEN BEGIN
FOR j := l to n-1 DO BEGIN
s := 0.0;
FOR k := i to m-1 DO BEGIN
s := s+a[k,i]*a[k,j]
END;
f := s/h;
FOR k := i to m-1 DO BEGIN
a[k,j] := a[k,j]+
f*a[k,i]
END
END
END;
FOR k := i to m-1 DO BEGIN
a[k,i] := scale*a[k,i]
END
END
END;
w[i,i] := scale*g;
g := 0.0;
s := 0.0;
scale := 0.0;
IF ((i <= m-1) AND (i <> n-1)) THEN BEGIN
FOR k := l to n-1 DO BEGIN
scale := scale+abs(a[i,k])
END;
IF (scale <> 0.0) THEN BEGIN
FOR k := l to n-1 DO BEGIN
a[i,k] := a[i,k]/scale;
s := s+a[i,k]*a[i,k]
END;
f := a[i,l];
g := -sign(sqrt(s),f);
h := f*g-s;
a[i,l] := f-g;
FOR k := l to n-1 DO BEGIN
rv1[k] := a[i,k]/h
END;
IF (i <> m-1) THEN BEGIN
FOR j := l to m-1 DO BEGIN
s := 0.0;
FOR k := l to n-1 DO BEGIN
s := s+a[j,k]*a[i,k]
END;
FOR k := l to n-1 DO BEGIN
a[j,k] := a[j,k]
+s*rv1[k]
END
END
END;
FOR k := l to n-1 DO BEGIN
a[i,k] := scale*a[i,k]
END
END
END;
anorm := max(anorm,(abs(w[i,i])+abs(rv1[i])))
END;
FOR i := n-1 DOWNTO 0 DO BEGIN
IF (i < n-1) THEN BEGIN
IF (g <> 0.0) THEN BEGIN
FOR j := l to n-1 DO BEGIN
v[j,i] := (a[i,j]/a[i,l])/g
END;
FOR j := l to n-1 DO BEGIN
s := 0.0;
FOR k := l to n-1 DO BEGIN
s := s+a[i,k]*v[k,j]
END;
FOR k := l to n-1 DO BEGIN
v[k,j] := v[k,j]+s*v[k,i]
END
END
END;
FOR j := l to n-1 DO BEGIN
v[i,j] := 0.0;
v[j,i] := 0.0
END
END;
v[i,i] := 1.0;
g := rv1[i];
l := i
END;
FOR i := n-1 DOWNTO 0 DO BEGIN
l := i+1;
g := w[i,i];
IF (i < n-1) THEN BEGIN
FOR j := l to n-1 DO BEGIN
a[i,j] := 0.0
END
END;
IF (g <> 0.0) THEN BEGIN
g := 1.0/g;
IF (i <> n-1) THEN BEGIN
FOR j := l to n-1 DO BEGIN
s := 0.0;
FOR k := l to m-1 DO BEGIN
s := s+a[k,i]*a[k,j]
END;
f := (s/a[i,i])*g;
FOR k := i to m-1 DO BEGIN
a[k,j] := a[k,j]+f*a[k,i]
END
END
END;
FOR j := i to m-1 DO BEGIN
a[j,i] := a[j,i]*g
END
END ELSE BEGIN
FOR j := i to m-1 DO BEGIN
a[j,i] := 0.0
END
END;
a[i,i] := a[i,i]+1.0
END;
FOR k := n-1 DOWNTO 0 DO BEGIN
FOR its := 1 to 30 DO BEGIN
FOR l := k DOWNTO 0 DO BEGIN
nm := l-1;
IF ((abs(rv1[l])+anorm) = anorm) THEN GOTO 2;
IF ((abs(w[nm,nm])+anorm) = anorm) THEN GOTO 1
END;
1: c := 0.0;
s := 1.0;
FOR i := l to k DO BEGIN
f := s*rv1[i];
IF ((abs(f)+anorm) <> anorm) THEN BEGIN
g := w[i,i];
h := sqrt(f*f+g*g);
w[i,i] := h;
h := 1.0/h;
c := (g*h);
s := -(f*h);
FOR j := 0 to m-1 DO BEGIN
y := a[j,nm];
z := a[j,i];
a[j,nm] := (y*c)+(z*s);
a[j,i] := -(y*s)+(z*c)
END
END
END;
2: z := w[k,k];
IF (l = k) THEN BEGIN
IF (z < 0.0) THEN BEGIN
w[k,k] := -z;
FOR j := 0 to n-1 DO BEGIN
v[j,k] := -v[j,k]
END
END;
GOTO 3
END;
IF (its = 30) THEN BEGIN
showmessage('No convergence in 30 SVDCMP iterations');
exit;
END;
x := w[l,l];
nm := k-1;
y := w[nm,nm];
g := rv1[nm];
h := rv1[k];
f := ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g := sqrt(f*f+1.0);
f := ((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
c := 1.0;
s := 1.0;
FOR j := l to nm DO BEGIN
i := j+1;
g := rv1[i];
y := w[i,i];
h := s*g;
g := c*g;
z := sqrt(f*f+h*h);
rv1[j] := z;
c := f/z;
s := h/z;
f := (x*c)+(g*s);
g := -(x*s)+(g*c);
h := y*s;
y := y*c;
FOR nm := 0 to n-1 DO BEGIN
x := v[nm,j];
z := v[nm,i];
v[nm,j] := (x*c)+(z*s);
v[nm,i] := -(x*s)+(z*c)
END;
z := sqrt(f*f+h*h);
w[j,j] := z;
IF (z <> 0.0) THEN BEGIN
z := 1.0/z;
c := f*z;
s := h*z
END;
f := (c*g)+(s*y);
x := -(s*g)+(c*y);
FOR nm := 0 to m-1 DO BEGIN
y := a[nm,j];
z := a[nm,i];
a[nm,j] := (y*c)+(z*s);
a[nm,i] := -(y*s)+(z*c)
END
END;
rv1[l] := 0.0;
rv1[k] := f;
w[k,k] := x
END;
3: END;
{ mat_print(m,a,'U matrix');
mat_print(n,v,'V matrix');
writeln(lst,'Diagonal values of W inverse matrix');
for i := 1 to n do
write(lst,1/w[i]:6:3);
writeln(lst); }
for i := 0 to n-1 do
for j := 0 to n-1 do
begin
if w[i,i] < 1.0e-6 then vtimesw[i,j] := 0
else vtimesw[i,j] := v[i,j] * (1.0 / w[j,j] );
end;
{ mat_print(n,vtimesw,'V matrix times w inverse '); }
for i := 0 to m-1 do
for j := 0 to n-1 do
begin
ainverse[i,j] := 0.0;
for k := 0 to m-1 do
begin
ainverse[i,j] := ainverse[i,j] + vtimesw[i,k] * a[j,k]
end;
end;
{ mat_print(n,ainverse,'Inverse Matrix'); }
for i := 0 to n-1 do
for j := 0 to n-1 do
a[i,j] := ainverse[i,j];
ainverse := nil;
rv1 := nil;
end;
end.
Reference in New Issue View Git Blame Copy Permalink