203 lines
7.5 KiB
Plaintext
203 lines
7.5 KiB
Plaintext
BEGIN # Multiple Regression - trnslation of the VB.NET sample but using the #
|
|
# "to reduced row echelon form" routine from the Reduced row echelon Task #
|
|
|
|
PROC require = ( BOOL condition, STRING message )VOID:
|
|
IF NOT condition THEN
|
|
print( ( message, newline ) );
|
|
stop
|
|
FI # requiree # ;
|
|
|
|
MODE MATRIX = STRUCT( REF[,]REAL data
|
|
, INT row count
|
|
, INT col count
|
|
);
|
|
|
|
PRIO NEWMATRIX = 1;
|
|
OP NEWMATRIX = ( INT m rows, INT m cols )MATRIX:
|
|
BEGIN
|
|
MATRIX result;
|
|
require( m rows > 0, "Need at least one row" );
|
|
row count OF result := m rows;
|
|
require( m cols > 0, "Need at least one column" );
|
|
col count OF result := m cols;
|
|
data OF result := HEAP[ 1 : m rows, 1 : m cols ]REAL;
|
|
FOR r TO m rows DO FOR c TO m cols DO ( data OF result )[ r, c ] := 0 OD OD;
|
|
result
|
|
END # NEWMATRIX # ;
|
|
|
|
OP NEWMATRIX = ( [,]REAL source )MATRIX:
|
|
BEGIN
|
|
MATRIX result;
|
|
INT m rows = 1 + ( 1 UPB source - 1 LWB source );
|
|
require( m rows > 0, "Need at least one row" );
|
|
row count OF result := m rows;
|
|
INT m cols = 1 + ( 2 UPB source - 2 LWB source );
|
|
require( m cols > 0, "Need at least one column" );
|
|
col count OF result := m cols;
|
|
data OF result := HEAP[ 1 : m rows, 1 : m cols ]REAL := source[ AT 1, AT 1 ];
|
|
result
|
|
END # NEWMATRIX # ;
|
|
|
|
OP NEWMATRIX = ( []REAL source )MATRIX: # New Matrix(ConvertArray(source)) #
|
|
BEGIN
|
|
INT len = 1 + ( UPB source - LWB source );
|
|
[ 1 : 1, 1 : len ]REAL dest;
|
|
dest[ 1, : ] := source;
|
|
NEWMATRIX dest
|
|
END # NEWMATRIX # ;
|
|
|
|
OP * = ( MATRIX m1, m2 )MATRIX:
|
|
BEGIN
|
|
INT rc1 = row count OF m1;
|
|
INT cc1 = col count OF m1;
|
|
INT rc2 = row count OF m2;
|
|
INT cc2 = col count OF m2;
|
|
require( cc1 = rc2, "Cannot multiply if the first columns does not equal the second rows" );
|
|
MATRIX result := rc1 NEWMATRIX cc2;
|
|
FOR i TO rc1 DO
|
|
FOR j TO cc2 DO
|
|
FOR k TO rc2 DO
|
|
( data OF result ) [ i, j ] +:= ( data OF m1 )[ i, k ]
|
|
* ( data OF m2 )[ k, j ]
|
|
OD
|
|
OD
|
|
OD;
|
|
result
|
|
END # * # ;
|
|
|
|
PROC transpose = ( MATRIX m )MATRIX:
|
|
BEGIN
|
|
INT rc = row count OF m;
|
|
INT cc = col count OF m;
|
|
MATRIX trans := cc NEWMATRIX rc;
|
|
FOR i TO cc DO
|
|
FOR j TO rc DO
|
|
( data OF trans )[ i, j ] := ( data OF m )[ j, i ]
|
|
OD
|
|
OD;
|
|
trans
|
|
END # transpose # ;
|
|
|
|
# BEGIN code from the Reduced row echelon form task #
|
|
MODE FIELD = REAL; # FIELD can be REAL, LONG REAL etc, or COMPL, FRAC etc #
|
|
MODE VEC = [0]FIELD;
|
|
MODE MAT = [0,0]FIELD;
|
|
PROC to reduced row echelon form = (REF MAT m)VOID: (
|
|
INT lead col := 2 LWB m;
|
|
|
|
FOR this row FROM LWB m TO UPB m DO
|
|
IF lead col > 2 UPB m THEN return FI;
|
|
INT other row := this row;
|
|
WHILE m[other row,lead col] = 0 DO
|
|
other row +:= 1;
|
|
IF other row > UPB m THEN
|
|
other row := this row;
|
|
lead col +:= 1;
|
|
IF lead col > 2 UPB m THEN return FI
|
|
FI
|
|
OD;
|
|
IF this row /= other row THEN
|
|
VEC swap = m[this row,lead col:];
|
|
m[this row,lead col:] := m[other row,lead col:];
|
|
m[other row,lead col:] := swap
|
|
FI;
|
|
FIELD scalef = 1/m[this row,lead col];
|
|
IF scalef /= 1 THEN
|
|
m[this row,lead col] := 1;
|
|
FOR col FROM lead col+1 TO 2 UPB m DO m[this row,col] *:= scalef OD
|
|
FI;
|
|
FOR other row pos FROM LWB m TO UPB m DO
|
|
IF this row /= other row pos THEN
|
|
REAL o scale = m[other row pos,lead col];
|
|
m[other row pos,lead col]:=0;
|
|
FOR col FROM lead col+1 TO 2 UPB m DO m[other row pos,col] -:= o scale*m[this row,col] OD
|
|
FI
|
|
OD;
|
|
lead col +:= 1
|
|
OD;
|
|
return: EMPTY
|
|
);
|
|
# END code from the Reduced row echelon form task #
|
|
|
|
PROC inverse = ( MATRIX m )MATRIX:
|
|
BEGIN
|
|
require( row count OF m = col count OF m, "Not a square matrix" );
|
|
INT len = row count OF m;
|
|
MATRIX aug := len NEWMATRIX ( 2 * len );
|
|
FOR i TO len DO
|
|
FOR j TO len DO
|
|
( data OF aug )[ i, j ] := ( data OF m )[ i, j ];
|
|
( data OF aug )[ i, j + len ] := 0
|
|
OD;
|
|
# augment identity matrix to right #
|
|
( data OF aug )[ i, i + len ] := 1.0
|
|
OD;
|
|
to reduced row echelon form( data OF aug );
|
|
MATRIX inv := len NEWMATRIX len;
|
|
FOR i TO len DO
|
|
FOR j FROM len + 1 TO 2 * len DO
|
|
( data OF inv)[ i, j - len ] := ( data OF aug )[ i, j ]
|
|
OD
|
|
OD;
|
|
inv
|
|
END # inverse # ;
|
|
|
|
PROC multiple regression = ( []REAL y, MATRIX x )[]REAL:
|
|
BEGIN
|
|
MATRIX tm := NEWMATRIX y;
|
|
MATRIX cy := NEWMATRIX data OF transpose( tm );
|
|
MATRIX cx := NEWMATRIX data OF transpose( x );
|
|
( data OF transpose( inverse( x * cx ) * x * cy ) )[ 1, : ]
|
|
END # multiple regression # ;
|
|
|
|
OP PRINTARRAY = ( []REAL list )VOID:
|
|
BEGIN
|
|
print( ( "[" ) );
|
|
FOR i FROM LWB list TO UPB list DO
|
|
# convert list[ i ] to a string, remove trailing 0s and leading spaces #
|
|
STRING v := fixed( list[ i ], -20, 15 )[ AT 1 ];
|
|
WHILE v[ UPB v ] = "0" DO v := v[ : UPB v - 1 ] OD;
|
|
IF v[ UPB v ] = "." THEN v := v[ : UPB v - 1 ] FI;
|
|
WHILE v[ 1 ] = " " DO v := v[ 2 : ] OD;
|
|
print( ( IF i > LWB list THEN ", " ELSE "" FI, v ) )
|
|
OD;
|
|
print( ( "]" ) )
|
|
END # PRINTARRAY # ;
|
|
|
|
BEGIN
|
|
[]REAL y = ( 1.0, 2.0, 3.0, 4.0, 5.0 );
|
|
MATRIX x := NEWMATRIX []REAL( 2.0, 1.0, 3.0, 4.0, 5.0 );
|
|
[]REAL v = multiple regression( y, x );
|
|
PRINTARRAY v;
|
|
print( ( newline ) )
|
|
END;
|
|
BEGIN
|
|
[]REAL y = ( 3.0, 4.0, 5.0 );
|
|
MATRIX x := NEWMATRIX [,]REAL( ( 1.0, 2.0, 1.0 )
|
|
, ( 1.0, 1.0, 2.0 )
|
|
);
|
|
[]REAL v = multiple regression( y, x );
|
|
PRINTARRAY v;
|
|
print( ( newline ) )
|
|
END;
|
|
BEGIN
|
|
[]REAL y = ( 52.21, 53.12, 54.48, 55.84, 57.2, 58.57, 59.93
|
|
, 61.29, 63.11, 64.47, 66.28, 68.1, 69.92, 72.19, 74.46
|
|
);
|
|
[]REAL a = ( 1.47, 1.5, 1.52, 1.55, 1.57, 1.6, 1.63, 1.65
|
|
, 1.68, 1.7, 1.73, 1.75, 1.78, 1.8, 1.83
|
|
);
|
|
[ 1 : 3, 1 : 1 + ( UPB a - LWB a ) ]REAL xs;
|
|
FOR i FROM LWB a TO UPB a DO
|
|
xs[ 1, i ] := 1.0;
|
|
xs[ 2, i ] := a[ i ];
|
|
xs[ 3, i ] := a[ i ] * a[ i ]
|
|
OD;
|
|
MATRIX x := NEWMATRIX xs;
|
|
[]REAL v = multiple regression( y, x );
|
|
PRINTARRAY v;
|
|
print( ( newline ) )
|
|
END
|
|
|
|
END
|