Multiple regression: Difference between revisions
Added various BASIC dialects (BASIC256, Chipmunk Basic, QBasic, QB64, True BASIC and Yabasic)
m (→{{header|R}}) |
(Added various BASIC dialects (BASIC256, Chipmunk Basic, QBasic, QB64, True BASIC and Yabasic)) |
||
(13 intermediate revisions by 8 users not shown) | |||
Line 20:
matrices.ads:
<
type Element_Type is private;
Zero : Element_Type;
Line 57:
Not_Square_Matrix : exception;
Not_Invertible : exception;
end Matrices;</
matrices.adb:
<
function "*" (Left, Right : Matrix) return Matrix is
Result : Matrix (Left'Range (1), Right'Range (2)) :=
Line 248:
return Result;
end Transpose;
end Matrices;</
Example multiple_regression.adb:
<
with Matrices;
procedure Multiple_Regression is
Line 318:
Output_Matrix (Float_Matrices.To_Matrix (Coefficients));
end;
end Multiple_Regression;</
{{out}}
Line 358:
-1.51161E+02
6.43514E+01</pre>
=={{header|ALGOL 68}}==
{{Trans|Visual Basic .NET}}...but using the "to reduced echelon form" routine from [[Reduced row echelon form#ALGOL_68]].
Algol 68 doesn't have classes, though it does have structures.
<br>Other than that, the main differences between this and the VB.NET sample are that Algol 68 has array slicing built in and generally uses a lower bound of 1 rather than 0 for arrays.
<br>Also, although <code>( ( 1, 2, 3 ), ( 6, 5, 4 ) )</code> is a 2x3 array, <code>( ( 1, 2, 3 ) )</code> is a 3x1 array (because <code>( x )</code> is not an array, so the outer brackets are superfluous, leaving <code>( 1, 2, 3 )</code> - i.e. a 1-dimensoional array - as the context requires a two-dimensional array, each value is coerced to an array resulting in a 3x1 two-dimensional array).
<syntaxhighlight lang="algol68">
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 rows, INT cols )MATRIX:
BEGIN
MATRIX result;
require( rows > 0, "Need at least one row" );
row count OF result := rows;
require( cols > 0, "Need at least one column" );
col count OF result := cols;
data OF result := HEAP[ 1 : rows, 1 : cols ]REAL;
FOR r TO rows DO FOR c TO cols DO ( data OF result )[ r, c ] := 0 OD OD;
result
END # NEWMATRIX # ;
OP NEWMATRIX = ( [,]REAL source )MATRIX:
BEGIN
MATRIX result;
INT rows = 1 + ( 1 UPB source - 1 LWB source );
require( rows > 0, "Need at least one row" );
row count OF result := rows;
INT cols = 1 + ( 2 UPB source - 2 LWB source );
require( cols > 0, "Need at least one column" );
col count OF result := cols;
data OF result := HEAP[ 1 : rows, 1 : 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 scale = 1/m[this row,lead col];
IF scale /= 1 THEN
m[this row,lead col] := 1;
FOR col FROM lead col+1 TO 2 UPB m DO m[this row,col] *:= scale OD
FI;
FOR other row FROM LWB m TO UPB m DO
IF this row /= other row THEN
REAL scale = m[other row,lead col];
m[other row,lead col]:=0;
FOR col FROM lead col+1 TO 2 UPB m DO m[other row,col] -:= 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
</syntaxhighlight>
{{out}}
<pre>
[0.981818181818182]
[1, 2]
[128.812803581374, -143.162022866676, 61.9603254433538]
</pre>
=={{header|BASIC}}==
==={{header|BASIC256}}===
{{trans|QB64}}
<syntaxhighlight lang="vbnet">arraybase 0
N = 14: M = 2: Q = 3 # number of points and M.R. polynom degree
dim X(N+1) # data points
X[0] = 1.47 : X[1] = 1.50 : X[2] = 1.52
X[3] = 1.55 : X[4] = 1.57 : X[5] = 1.60
X[6] = 1.63 : X[7] = 1.65 : X[8] = 1.68
X[9] = 1.70 : X[10] = 1.73 : X[11] = 1.75
X[12] = 1.78 : X[13] = 1.80 : X[14] = 1.83
dim Y(N+1) # data points
Y[0] = 52.21 : Y[1] = 53.12 : Y[2] = 54.48
Y[3] = 55.84 : Y[4] = 57.20 : Y[5] = 58.57
Y[6] = 59.93 : Y[7] = 61.29 : Y[8] = 63.11
Y[9] = 64.47 : Y[10] = 66.28 : Y[11] = 68.10
Y[12] = 69.92 : Y[13] = 72.19 : Y[14] = 74.46
dim S(N+1): dim T(N+1) # linear system coefficient
dim A(M+1, Q+1) # system to be solved
dim p(M+1, Q+1)
for k = 0 to 2 * M
S[k] = 0: T[k] = 0
for i = 0 to N
S[k] = S[k] + X[i] ^ k
if k <= M then T[k] = T[k] + Y[i] * X[i] ^ k
next i
next k
# build linear system
for fila = 0 to M
for columna = 0 to M
A[fila, columna] = S[fila + columna]
next columna
A[fila, columna] = T[fila]
next fila
print "Linear system coefficents:"
for i = 0 to M
for j = 0 to M + 1
print rjust(A[i, j], 14.1);
next j
print
next i
for j = 0 to M
for i = j to M
if A[i, j] <> 0 then exit for
next i
if i = M + 1 then print: print "SINGULAR MATRIX '" : end
for k = 0 to M + 1
p[j,k] = A[i,k] : A[i,k] = p[j,k] : A[j,k] = A[i,k]
next k
z = 1 / A[j, j]
for k = 0 to M + 1
A[j, k] = z * A[j, k]
next k
for i = 0 to M
if i <> j then
z = -A[i, j]
for k = 0 to M + 1
A[i, k] = A[i, k] + z * A[j, k]
next k
end if
next i
next j
print: print "Solutions:"
for i = 0 to M
print rjust(A[i, M + 1], 16.7);
next i
end</syntaxhighlight>
{{out}}
<pre>Linear system coefficents:
15.0 24.76 41.0532 931.17
24.76 41.0532 68.367934 1548.2453
41.0532 68.367934 114.34828728 2585.543151
Solutions:
128.81280358 -143.162022867 61.9603254432</pre>
==={{header|Chipmunk Basic}}===
{{trans|FreeBASIC}}
{{works with|Chipmunk Basic|3.6.4}}
<syntaxhighlight lang="vbnet">100 cls
110 n = 14 : m = 2 : q = 3 ' number of points and M.R. polynom degree
120 dim x(n)' data points
130 data 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
140 for c = 0 to n
150 read x(c)
160 next c
170 dim y(n)' data points
180 data 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
190 for c = 0 to n
200 read y(c)
210 next c
220 dim s(n) : dim t(n)' linear system coefficient
230 dim a(m,q)' sistem to be solved
240 dim p(m,q)
250 for k = 0 to 2*m
260 s(k) = 0 : t(k) = 0
270 for i = 0 to n
280 s(k) = s(k)+x(i)^k
290 if k <= m then t(k) = t(k)+y(i)*x(i)^k
300 next i
310 next k
320 ' build linear system
330 for fila = 0 to m
340 for columna = 0 to m
350 a(fila,columna) = s(fila+columna)
360 next columna
370 a(fila,columna) = t(fila)
380 next fila
390 print "Linear system coefficents:"
400 for i = 0 to m
410 for j = 0 to m+1
420 print using "######.#";a(i,j);
430 next j
440 print
450 next i
460 for j = 0 to m
470 for i = j to m
480 if a(i,j) <> 0 then goto 500
490 next i
500 if i = m+1 then
510 print : print "SINGULAR MATRIX '"
520 end
530 endif
540 for k = 0 to m+1
560 p(j,k) = a(i,k)
570 a(i,k) = p(j,k)
580 a(j,k) = a(i,k)
590 next k
600 z = 1/a(j,j)
610 for k = 0 to m+1
620 a(j,k) = z*a(j,k)
630 next k
640 for i = 0 to m
650 if i <> j then
660 z = -a(i,j)
670 for k = 0 to m+1
680 a(i,k) = a(i,k)+z*a(j,k)
690 next k
700 endif
710 next i
720 next j
730 print : print "Solutions:"
740 for i = 0 to m
750 print using " #####.#######";a(i,m+1);
760 next i
770 print
780 end</syntaxhighlight>
{{out}}
<pre>Same as FreeBASIC entry.</pre>
==={{header|QBasic}}===
{{works with|QBasic|1.1}}
{{works with|QuickBasic|4.5}}
The [[#QB64|QB64]] solution works without any changes.
==={{header|QB64}}===
{{trans|FreeBASIC}}
<syntaxhighlight lang="vbnet">Const N = 14, M = 2, Q = 3 ' number of points and M.R. polynom degree
Dim X(N) As Double ' data points
Data 1.47,1.50,1.52,1.55,1.57,1.60,1.63,1.65,1.68,1.70,1.73,1.75,1.78,1.80,1.83
For c = LBound(X) To UBound(X)
Read X(c)
Next c
Dim As Double Y(N) ' data points
Data 52.21,53.12,54.48,55.84,57.20,58.57,59.93,61.29,63.11,64.47,66.28,68.10,69.92,72.19,74.46
For c = LBound(Y) To UBound(Y)
Read Y(c)
Next c
Dim As Double S(N), T(N) ' linear system coefficient
Dim As Double A(M, Q) ' system to be solved
Dim As Integer i, k, j, fila, columna
Dim As Double z
For k = 0 To 2 * M
S(k) = 0: T(k) = 0
For i = 0 To N
S(k) = S(k) + X(i) ^ k
If k <= M Then T(k) = T(k) + Y(i) * X(i) ^ k
Next i
Next k
' build linear system
For fila = 0 To M
For columna = 0 To M
A(fila, columna) = S(fila + columna)
Next columna
A(fila, columna) = T(fila)
Next fila
Print "Linear system coefficents:"
For i = 0 To M
For j = 0 To M + 1
Print Using "######.#"; A(i, j);
Next j
Print
Next i
For j = 0 To M
For i = j To M
If A(i, j) <> 0 Then Exit For
Next i
If i = M + 1 Then
Print: Print "SINGULAR MATRIX '"
End
End If
For k = 0 To M + 1
Swap A(j, k), A(i, k)
Next k
z = 1 / A(j, j)
For k = 0 To M + 1
A(j, k) = z * A(j, k)
Next k
For i = 0 To M
If i <> j Then
z = -A(i, j)
For k = 0 To M + 1
A(i, k) = A(i, k) + z * A(j, k)
Next k
End If
Next i
Next j
Print: Print "Solutions:"
For i = 0 To M
Print Using " #####.#######"; A(i, M + 1);
Next i
End</syntaxhighlight>
{{out}}
<pre>Same as FreeBASIC entry.</pre>
==={{header|True BASIC}}===
{{trans|QB64}}
<syntaxhighlight lang="qbasic">OPTION BASE 0
LET n = 14 ! number of points and M.R. polynom degree
LET m = 2
LET q = 3
DIM x(0) ! data points
MAT REDIM x(n)
DATA 1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83
FOR c = LBOUND(x) TO UBOUND(x)
READ x(c)
NEXT c
DIM y(0) ! data points
MAT REDIM y(n)
DATA 52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46
FOR c = LBOUND(y) TO UBOUND(y)
READ y(c)
NEXT c
DIM s(0) ! linear system coefficient
MAT REDIM s(n)
DIM t(0)
MAT REDIM t(n)
DIM a(0,0) ! system to be solved
MAT REDIM a(m, q)
DIM p(0,0)
MAT REDIM p(m, q)
FOR k = 0 TO 2*m
LET s(k) = 0
LET t(k) = 0
FOR i = 0 TO n
LET s(k) = s(k)+x(i)^k
IF k <= m THEN LET t(k) = t(k)+y(i)*x(i)^k
NEXT i
NEXT k
! build linear system
FOR fila = 0 TO m
FOR columna = 0 TO m
LET a(fila, columna) = s(fila+columna)
NEXT columna
LET a(fila, columna) = t(fila)
NEXT fila
PRINT "Linear system coefficents:"
FOR i = 0 TO m
FOR j = 0 TO m+1
PRINT USING "######.#": a(i, j);
NEXT j
PRINT
NEXT i
FOR j = 0 TO m
FOR i = j TO m
IF a(i, j) <> 0 THEN EXIT FOR
NEXT i
IF i = m+1 THEN
PRINT
PRINT "SINGULAR MATRIX '"
STOP
END IF
FOR k = 0 TO m+1
LET p(j, k) = a(i, k)
LET a(i, k) = p(j, k)
LET a(j, k) = a(i, k)
NEXT k
LET z = 1/a(j, j)
FOR k = 0 TO m+1
LET a(j, k) = z*a(j, k)
NEXT k
FOR i = 0 TO m
IF i <> j THEN
LET z = -a(i, j)
FOR k = 0 TO m+1
LET a(i, k) = a(i, k)+z*a(j, k)
NEXT k
END IF
NEXT i
NEXT j
PRINT
PRINT "Solutions:"
FOR i = 0 TO m
PRINT USING " #####.#######": a(i, m+1);
NEXT i
END</syntaxhighlight>
==={{header|Yabasic}}===
{{trans|QB64}}
<syntaxhighlight lang="vbnet">// number of points and M.R. polynom degree
N = 14
M = 2
Q = 3
dim X(N) // data points
X(0) = 1.47 : X(1) = 1.50 : X(2) = 1.52
X(3) = 1.55 : X(4) = 1.57 : X(5) = 1.60
X(6) = 1.63 : X(7) = 1.65 : X(8) = 1.68
X(9) = 1.70 : X(10) = 1.73 : X(11) = 1.75
X(12) = 1.78 : X(13) = 1.80 : X(14) = 1.83
dim Y(N) // data points
Y(0) = 52.21 : Y(1) = 53.12 : Y(2) = 54.48
Y(3) = 55.84 : Y(4) = 57.20 : Y(5) = 58.57
Y(6) = 59.93 : Y(7) = 61.29 : Y(8) = 63.11
Y(9) = 64.47 : Y(10) = 66.28 : Y(11) = 68.10
Y(12) = 69.92 : Y(13) = 72.19 : Y(14) = 74.46
dim S(N), T(N) // linear system coefficient
dim A(M, Q) // sistem to be solved
dim p(M, Q)
for k = 0 to 2*M
S(k) = 0 : T(k) = 0
for i = 0 to N
S(k) = S(k) + X(i) ^ k
if k <= M T(k) = T(k) + Y(i) * X(i) ^ k
next i
next k
// build linear system
for fila = 0 to M
for columna = 0 to M
A(fila, columna) = S(fila+columna)
next columna
A(fila, columna) = T(fila)
next fila
print "Linear system coefficents:"
for i = 0 to M
for j = 0 to M+1
print A(i,j) using "#####.#";
next j
print
next i
for j = 0 to M
for i = j to M
if A(i,j) <> 0 break
next i
if i = M+1 then
print "\nSINGULAR MATRIX '"
end
end if
for k = 0 to M+1
p(j,k) = A(i,k) : A(i,k) = p(j,k) : A(j,k) = A(i,k)
next k
z = 1 / A(j,j)
for k = 0 to M+1
A(j,k) = z * A(j,k)
next k
for i = 0 to M
if i <> j then
z = -A(i,j)
for k = 0 to M+1
A(i,k) = A(i,k) + z * A(j,k)
next k
end if
next i
next j
print "\nSolutions:"
for i = 0 to M
print A(i,M+1) using "#######.#######";
next i
print
end</syntaxhighlight>
=={{header|BBC BASIC}}==
{{works with|BBC BASIC for Windows}}
<
INSTALL @lib$+"ARRAYLIB"
Line 389 ⟶ 1,005:
t() = t().m()
c() = y().t()
ENDPROC</
{{out}}
<pre>
Line 396 ⟶ 1,012:
=={{header|C}}==
Using GNU gsl and c99, with the WP data<
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
Line 439 ⟶ 1,055:
gsl_multifit_linear_free(wspc);
}</
=={{header|C++}}==
{{trans|Java}}
<
#include <iostream>
Line 669 ⟶ 1,285:
return 0;
}</
{{out}}
<pre>[0.981818]
Line 677 ⟶ 1,293:
=={{header|C sharp|C#}}==
{{libheader|Math.Net}}
<
using MathNet.Numerics.LinearRegression;
using MathNet.Numerics.LinearAlgebra;
Line 694 ⟶ 1,310:
Console.WriteLine(β);
}
}</
{{out}}
Line 706 ⟶ 1,322:
Uses the routine (chol A) from [[Cholesky decomposition]], (mmul A B) from [[Matrix multiplication]], (mtp A) from [[Matrix transposition]].
<
;; Solve a linear system AX=B where A is symmetric and positive definite, so it can be Cholesky decomposed.
(defun linsys (A B)
Line 740 ⟶ 1,356:
(linsys (mmul (mtp A) A)
(mmul (mtp A) b)))
</syntaxhighlight>
To show an example of multiple regression, (polyfit x y n) from [[Polynomial regression]], which itself uses (linsys A B) and (lsqr A b), will be used to fit a second degree order polynomial to data.
<
(y (make-array '(1 11) :initial-contents '((1 6 17 34 57 86 121 162 209 262 321)))))
(polyfit x y 2))
#2A((0.9999999999999759d0) (2.000000000000005d0) (3.0d0))</
=={{header|D}}==
{{trans|Java}}
<
import std.array;
import std.exception;
Line 975 ⟶ 1,591:
v = multipleRegression(y, x);
v.writeln;
}</
{{out}}
<pre>[0.981818]
Line 985 ⟶ 1,601:
{{libheader|calc}}
<
(x2 '(0 1 1 3 3 7 6 7 3 9 8))
(y '(1 6 17 34 57 86 121 162 209 262 321)))
(apply #'calc-eval "fit(a*X1+b*X2+c,[X1,X2],[a,b,c],[$1 $2 $3])" nil
(mapcar (lambda (items) (cons 'vec items)) (list x1 x2 y))))</
{{out}}
Line 996 ⟶ 1,612:
=={{header|ERRE}}==
<
!$DOUBLE
Line 1,077 ⟶ 1,693:
END FOR
END PROGRAM</
{{out}}
<pre>LINEAR SYSTEM COEFFICENTS
Line 1,097 ⟶ 1,713:
{{libheader|SLATEC}} [http://netlib.org/slatec/ Available at the Netlib]
<
* MR - multiple regression using the SLATEC library routine DHFTI
*
Line 1,182 ⟶ 1,798:
STOP 'program complete'
END
</syntaxhighlight>
{{out}}
<pre>
Line 1,188 ⟶ 1,804:
STOP program complete
</pre>
=={{header|FreeBASIC}}==
{{trans|ERRE}}
<syntaxhighlight lang="vb">Const N = 14, M = 2, Q = 3 ' number of points and M.R. polynom degree
Dim As Double X(0 to N) = {1.47,1.50,1.52,1.55,1.57, _
1.60,1.63,1.65,1.68,1.70,1.73,1.75,1.78,1.80,1.83} ' data points
Dim As Double Y(0 to N) = {52.21,53.12,54.48,55.84,57.20, _
58.57,59.93,61.29,63.11,64.47,66.28,68.10,69.92,72.19,74.46} ' data points
Dim As Double S(N), T(N) ' linear system coefficient
Dim As Double A(M, Q) ' sistem to be solved
Dim As Integer i, k, j, fila, columna
Dim as Double z
For k = 0 To 2*M
S(k) = 0 : T(k) = 0
For i = 0 To N
S(k) += X(i) ^ k
If k <= M Then T(k) += Y(i) * X(i) ^ k
Next i
Next k
' build linear system
For fila = 0 To M
For columna = 0 To M
A(fila, columna) = S(fila+columna)
Next columna
A(fila, columna) = T(fila)
Next fila
Print "Linear system coefficents:"
For i = 0 To M
For j = 0 To M+1
Print Using "######.#"; A(i,j);
Next j
Print
Next i
For j = 0 To M
For i = j To M
If A(i,j) <> 0 Then Exit For
Next i
If i = M+1 Then
Print !"\nSINGULAR MATRIX '"
Sleep: End
End If
For k = 0 To M+1
Swap A(j,k), A(i,k)
Next k
z = 1 / A(j,j)
For k = 0 To M+1
A(j,k) = z * A(j,k)
Next k
For i = 0 To M
If i <> j Then
z = -A(i,j)
For k = 0 To M+1
A(i,k) += z * A(j,k)
Next k
End If
Next i
Next j
Print !"\nSolutions:"
For i = 0 To M
Print Using " #####.#######"; A(i,M+1);
Next i
Sleep</syntaxhighlight>
{{out}}
<pre>Linear system coefficents:
15.0 24.8 41.1 931.2
24.8 41.1 68.4 1548.2
41.1 68.4 114.3 2585.5
Solutions:
128.8128036 -143.1620229 61.9603254</pre>
=={{header|Go}}==
The [http://en.wikipedia.org/wiki/Ordinary_least_squares#Example_with_real_data example] on WP happens to be a polynomial regression example, and so code from the [[Polynomial regression]] task can be reused here. The only difference here is that givens x and y are computed in a separate function as a task prerequisite.
===Library gonum/matrix===
<
import (
Line 1,224 ⟶ 1,917:
x, y := givens()
fmt.Printf("%.4f\n", mat64.Formatted(mat64.QR(x).Solve(y)))
}</
{{out}}
<pre>
Line 1,232 ⟶ 1,925:
</pre>
===Library go.matrix===
<
import (
Line 1,277 ⟶ 1,970:
}
fmt.Println(c)
}</
{{out}}
<pre>
Line 1,285 ⟶ 1,978:
=={{header|Haskell}}==
Using package [http://hackage.haskell.org/package/hmatrix hmatrix] from HackageDB
<
import Numeric.LinearAlgebra.LAPACK
Line 1,296 ⟶ 1,989:
v :: Matrix Double
v = (3><1)
[1.745005,-4.448092,-4.160842]</
Using lapack::dgels
<
(3><1)
[ 0.9335611922087276
, 1.101323491272865
, 1.6117769115824 ]</
Or
<
(3><1)
[ 0.9335611922087278
, 1.101323491272865
, 1.6117769115824006 ]</
=={{header|Hy}}==
<
[numpy [ones column-stack]]
[numpy.random [randn]]
Line 1,323 ⟶ 2,016:
(print (first (lstsq
(column-stack (, (ones n) x1 x2 (* x1 x2)))
y)))</
=={{header|J}}==
<
x=: 1.47 1.50 1.52 1.55 1.57 1.60 1.63 1.65 1.68 1.70 1.73 1.75 1.78 1.80 1.83
y=: 52.21 53.12 54.48 55.84 57.20 58.57 59.93 61.29 63.11 64.47 66.28 68.10 69.92 72.19 74.46
y %. x ^/ i.3 NB. calculate coefficients b1, b2 and b3 for 2nd degree polynomial
128.813 _143.162 61.9603</
Breaking it down:
<
X=: (x^0) ,. (x^1) ,. (x^2) NB. equivalent of previous line
4{.X NB. show first 4 rows of X
Line 1,346 ⟶ 2,039:
NB. y %. X does matrix division and gives the regression coefficients
y %. X
128.813 _143.162 61.9603</
In other words <tt> beta=: y %. X </tt> is the equivalent of:<br>
<math> \hat\beta = (X'X)^{-1}X'y</math><br>
To confirm:
<
NB. %.X is matrix inverse of X
NB. |:X is transpose of X
Line 1,360 ⟶ 2,053:
X (%.@:xpy@[ mp xpy) y
128.814 _143.163 61.9606
</syntaxhighlight>
LAPACK routines are also available via the Addon <tt>math/lapack</tt>.
<
load 'math/lapack/gels'
gels_jlapack_ X;y
128.813 _143.162 61.9603</
=={{header|Java}}==
{{trans|Kotlin}}
<
import java.util.Objects;
Line 1,578 ⟶ 2,271:
printVector(v);
}
}</
{{out}}
<pre>[0.9818181818181818]
Line 1,593 ⟶ 2,286:
Extends the Matrix class from [[Matrix Transpose#JavaScript]], [[Matrix multiplication#JavaScript]], [[Reduced row echelon form#JavaScript]].
Uses the IdentityMatrix from [[Matrix exponentiation operator#JavaScript]]
<
Matrix.prototype.inverse = function() {
if (this.height != this.width) {
Line 1,639 ⟶ 2,332:
)
);
print(y.regression_coefficients(x));</
{{out}}
<pre>0.9818181818181818
Line 1,656 ⟶ 2,349:
'''Preliminaries'''
<
reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );
Line 1,667 ⟶ 2,360:
reduce range(0; $p) as $j
(.;
.[$i][$j] = dot_product( A[$i]; $BT[$j] ) ));</
'''Multiple Regression'''
<
(y|transpose) as $cy
| (x|transpose) as $cx
Line 1,693 ⟶ 2,386:
range(0; ys|length) as $i
| multipleRegression(ys[$i]; xs[$i])</
{{out}}
<pre>
Line 1,707 ⟶ 2,400:
As in Matlab, the backslash or slash operator (depending on the matrix ordering) can be used for solving this problem, for example:
<
y = [52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46]
X = [x.^0 x.^1 x.^2];
b = X \ y</
{{out}}
<pre>
Line 1,721 ⟶ 2,414:
=={{header|Kotlin}}==
As neither the JDK nor the Kotlin Standard Library has matrix operations built-in, we re-use functions written for various other tasks.
<
typealias Vector = DoubleArray
Line 1,841 ⟶ 2,534:
v = multipleRegression(y, x)
printVector(v)
}</
{{out}}
Line 1,856 ⟶ 2,549:
First build a random dataset.
<
X:=<ArrayTools[RandomArray](n,4,distribution=normal)|Vector(n,1,datatype=float[8])>:
Y:=X.<4.2,-2.8,-1.4,3.1,1.75>+convert(ArrayTools[RandomArray](n,1,distribution=normal),Vector):</
Now the linear regression, with either the LinearAlgebra package, or the Statistics package.
<
# [4.33701132468683, -2.78654498997457, -1.41840666085642, 2.92065133466547, 1.76076442997642]
Line 1,880 ⟶ 2,573:
# R-squared: 0.9767, Adjusted R-squared: 0.9761
# 4.33701132468683 x1 - 2.78654498997457 x2 - 1.41840666085642 x3
# + 2.92065133466547 x4 + 1.76076442997642 c</
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
y = {52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46};
X = {x^0, x^1, x^2};
LeastSquares[Transpose@X, y]</
{{out}}
<pre>{128.813, -143.162, 61.9603}</pre>
Line 1,894 ⟶ 2,587:
The slash and backslash operator can be used for solving this problem. Here some random data is generated.
<
y = randn (1,n); % generate random vector y
X = randn (k,n); % generate random matrix X
b = y / X
b = 0.1457109 -0.0777564 -0.0712427 -0.0166193 0.0292955 -0.0079111 0.2265894 -0.0561589 -0.1752146 -0.2577663 </
In its transposed form yt = Xt * bt, the backslash operator can be used.
<
bt = Xt \ yt
bt =
Line 1,914 ⟶ 2,607:
-0.0561589
-0.1752146
-0.2577663</
Here is the example for estimating the polynomial fit
<
y = [52.21 53.12 54.48 55.84 57.20 58.57 59.93 61.29 63.11 64.47 66.28 68.10 69.92 72.19 74.46]
X = [x.^0;x.^1;x.^2];
b = y/X
128.813 -143.162 61.960</
Instead of "/", the slash operator, one can also write :
<
or
<
=={{header|Nim}}==
{{libheader|arraymancer}}
<
import math
Line 1,948 ⟶ 2,641:
var a = stack(height.ones_like, height, height *. height, axis = 1)
echo toSeq(least_squares_solver(a, weight).solution.items)</
{{out}}
Line 1,954 ⟶ 2,647:
=={{header|PARI/GP}}==
<
addhelp(pseudoinv, "pseudoinv(M): Moore pseudoinverse of the matrix M.");
y*pseudoinv(X)</
=={{header|Perl}}==
<
use warnings;
use Statistics::Regression;
Line 1,972 ⟶ 2,665:
my @coeff = $reg->theta();
printf "%-6s %8.3f\n", $model[$_], $coeff[$_] for 0..@model-1;</
{{out}}
<pre>const 128.813
Line 1,980 ⟶ 2,673:
=={{header|Phix}}==
{{trans|ERRE}}
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">N</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">15</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">M</span><span style="color: #0000FF;">=</span><span style="color: #000000;">3</span>
Line 2,036 ⟶ 2,729:
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Solutions:\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">columnize</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">M</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<!--</
{{out}}
<pre>
Line 2,048 ⟶ 2,741:
=={{header|PicoLisp}}==
<
# Matrix transposition
Line 2,090 ⟶ 2,783:
(car X) ) ) ) )
(T (> (inc 'Lead) Cols)) ) )
Mat )</
{{trans|JavaScript}}
<
(let N (length Mat)
(unless (= N (length (car Mat)))
Line 2,110 ⟶ 2,803:
X (columnVector (2.0 1.0 3.0 4.0 5.0)) )
(round (caar (regressionCoefficients Y X)) 17)</
{{out}}
<pre>-> "0.98181818181818182"</pre>
Line 2,117 ⟶ 2,810:
{{libheader|NumPy}}
'''Method with matrix operations'''
<
height = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,
Line 2,127 ⟶ 2,820:
y = np.mat(weight)
print(y * X.T * (X*X.T).I)</
{{out}}
<pre>
Line 2,133 ⟶ 2,826:
</pre>
'''Using numpy lstsq function'''
<
height = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,
Line 2,143 ⟶ 2,836:
y = weight
print(np.linalg.lstsq(X, y)[0])</
{{out}}
<pre>
Line 2,153 ⟶ 2,846:
R provides the '''lm''' function for linear regression.
<
y <- c(52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46)
lm( y ~ x + I(x^2))</
{{out}}
Line 2,170 ⟶ 2,863:
is useful to illustrate R's model description and linear algebra capabilities.
<
## parse and evaluate the model formula
Line 2,185 ⟶ 2,878:
}
simpleMultipleReg(y ~ x + I(x^2))</
This produces the same coefficients as lm()
Line 2,199 ⟶ 2,892:
than the method above, is to solve the linear system directly
and use the crossprod function:
<
A numerically more stable way is to use the QR decomposition of the design matrix:
<
mf <- model.frame(formula)
X <- model.matrix(mf)
Line 2,223 ⟶ 2,916:
# (Intercept) x I(x^2)
# 128.81280 -143.16202 61.96033</
=={{header|Racket}}==
<
#lang racket
(require math)
Line 2,233 ⟶ 2,926:
(define (fit X y)
(matrix-solve (matrix* (T X) X) (matrix* (T X) y)))
</syntaxhighlight>
Test:
<
(fit (matrix [[1 2]
[2 5]
Line 2,246 ⟶ 2,939:
{{out}}
(array #[#[9 1/3] #[-3 1/3]])
</syntaxhighlight>
=={{header|Raku}}==
{{broken}}
(formerly Perl 6)
We're going to solve the example on the Wikipedia article using [https://github.com/grondilu/clifford Clifford], a [https://en.wikipedia.org/wiki/Geometric_algebra geometric algebra] module. Optimization for large vector space does not quite work yet, so it's going to take (a lof of) time and a fair amount of memory, but it should work.
Line 2,276 ⟶ 2,970:
<syntaxhighlight lang="raku"
my @height = <1.47 1.50 1.52 1.55 1.57 1.60 1.63 1.65 1.68 1.70 1.73 1.75 1.78 1.80 1.83>;
my @weight = <52.21 53.12 54.48 55.84 57.20 58.57 59.93 61.29 63.11 64.47 66.28 68.10 69.92 72.19 74.46>;
Line 2,291 ⟶ 2,985:
say "α = ", ($w∧$h1∧$h2)·$I.reversion/$I2;
say "β = ", ($w∧$h2∧$h0)·$I.reversion/$I2;
say "γ = ", ($w∧$h0∧$h1)·$I.reversion/$I2;</
{{out}}
<pre>α = 128.81280357844
Line 2,304 ⟶ 2,998:
Using the standard library Matrix class:
<
def regression_coefficients y, x
Line 2,311 ⟶ 3,005:
(x.t * x).inverse * x.t * y
end</
Testing 2-dimension:
<
{{out}}
<pre>Matrix[[0.981818181818182]]</pre>
Line 2,320 ⟶ 3,014:
Testing 3-dimension:
Points(x,y,z): [1,1,3], [2,1,4] and [1,2,5]
<
{{out}}
<pre>Matrix[[0.9999999999999996], [2.0]]</pre>
Line 2,328 ⟶ 3,022:
First, build a random dataset:
<syntaxhighlight lang="text">set rng=mc seed=17760704.
new file.
input program.
Line 2,341 ⟶ 3,035:
end input program.
compute y=1.5+0.8*x1-0.7*x2+1.1*x3-1.7*x4+rv.normal(0,1).
execute.</
Now use the regression command:
<syntaxhighlight lang="text">regression /dependent=y
/method=enter x1 x2 x3 x4.</
{{out}}
<syntaxhighlight lang="text">Regression
Notes
|--------------------------------------------------------------------|---------------------------------------------------------------------------|
Line 2,426 ⟶ 3,120:
| |x4 |-1,770 |,073 |-,656 |-24,306|,000|
|----------------------------------------------------------------------------------------------|
a Dependent Variable: y</
=={{header|Stata}}==
Line 2,432 ⟶ 3,126:
First, build a random dataset:
<
set seed 17760704
set obs 200
Line 2,438 ⟶ 3,132:
gen x`i'=rnormal()
}
gen y=1.5+0.8*x1-0.7*x2+1.1*x3-1.7*x4+rnormal()</
Now, use the '''[https://www.stata.com/help.cgi?regress regress]''' command:
<syntaxhighlight lang
'''Output'''
Line 2,467 ⟶ 3,161:
The regress command also sets a number of '''[https://www.stata.com/help.cgi?ereturn ereturn]''' values, which can be used by subsequent commands. The coefficients and their standard errors also have a [https://www.stata.com/help.cgi?_variables special syntax]:
<
.75252466
Line 2,477 ⟶ 3,171:
. di _se[_cons]
.06978623</
See '''[https://www.stata.com/help.cgi?regress_postestimation regress postestimation]''' for a list of commands that can be used after a regression.
Line 2,483 ⟶ 3,177:
Here we compute [[wp:Akaike information criterion|Akaike's AIC]], the covariance matrix of the estimates, the predicted values and residuals:
<
Akaike's information criterion and Bayesian information criterion
Line 2,507 ⟶ 3,201:
. predict yhat, xb
. predict r, r</
=={{header|Tcl}}==
{{tcllib|math::linearalgebra}}
<
namespace eval multipleRegression {
namespace export regressionCoefficients
Line 2,526 ⟶ 3,220:
}
}
namespace import multipleRegression::regressionCoefficients</
Using an example from the Wikipedia page on the correlation of height and weight:
<
proc map {n exp list} {
upvar 1 $n v
Line 2,543 ⟶ 3,237:
}
# Wikipedia states that fitting up to the square of x[i] is worth it
puts [regressionCoefficients $y [map n {map v {expr {$v**$n}} $x} {0 1 2}]]</
{{out}} (a 3-vector of coefficients):
<pre>128.81280358170625 -143.16202286630732 61.96032544293041</pre>
Line 2,549 ⟶ 3,243:
=={{header|TI-83 BASIC}}==
{{works with|TI-83 BASIC|TI-84Plus 2.55MP}}
<
{52.21,53.12,54.48,55.84,57.20,58.57,59.93,61.29,63.11,64.47,66.28,68.10,69.92,72.19,74.46}→L₂
QuadReg L₁,L₂ </
{{out}}
<pre>
Line 2,564 ⟶ 3,258:
the Lapack library [http://www.netlib.org/lapack/lug/node27.html],
which is callable in Ursala like this:
<
test program:
<syntaxhighlight lang="ursala">x =
<
Line 2,577 ⟶ 3,271:
#cast %eL
example = regression_coefficients(x,y)</
The matrix x needn't be square, and has one row for each data point.
The length of y must equal the number of rows in x,
Line 2,594 ⟶ 3,288:
=={{header|Visual Basic .NET}}==
{{trans|Java}}
<
Sub Swap(Of T)(ByRef x As T, ByRef y As T)
Line 2,823 ⟶ 3,517:
End Sub
End Module</
{{out}}
<pre>[0.981818181818182]
Line 2,832 ⟶ 3,526:
{{trans|Kotlin}}
{{libheader|Wren-matrix}}
<
var multipleRegression = Fn.new { |y, x|
Line 2,859 ⟶ 3,553:
System.print(v)
System.print()
}</
{{out}}
Line 2,872 ⟶ 3,566:
=={{header|zkl}}==
Using the GNU Scientific Library:
<
height:=GSL.VectorFromData(1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,
1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83);
Line 2,880 ⟶ 3,574:
v.format().println();
GSL.Helpers.polyString(v).println();
GSL.Helpers.polyEval(v,height).format().println();</
{{out}}
<pre>
Line 2,890 ⟶ 3,584:
Or, using Lists:
{{trans|Common Lisp}}
<
fcn linsys(A,B){
n,m:=A.len(),B[1].len(); // A.rows,B.cols
Line 2,945 ⟶ 3,639:
if(M.len()==1) M[0].pump(List,List.create); // 1 row --> n columns
else M[0].zip(M.xplode(1));
}</
<
1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83));
weight:=T(T(52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93,
61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46));
polyfit(height,weight,2).flatten().println();</
{{out}}
<pre>
|