Multiple regression: Difference between revisions

m
mNo edit summary
m (→‎{{header|Wren}}: Minor tidy)
 
(38 intermediate revisions by 16 users not shown)
Line 20:
 
matrices.ads:
<langsyntaxhighlight Adalang="ada">generic
type Element_Type is private;
Zero : Element_Type;
Line 57:
Not_Square_Matrix : exception;
Not_Invertible : exception;
end Matrices;</langsyntaxhighlight>
 
matrices.adb:
<langsyntaxhighlight Adalang="ada">package body Matrices is
function "*" (Left, Right : Matrix) return Matrix is
Result : Matrix (Left'Range (1), Right'Range (2)) :=
Line 248:
return Result;
end Transpose;
end Matrices;</langsyntaxhighlight>
 
Example multiple_regression.adb:
<langsyntaxhighlight Adalang="ada">with Ada.Text_IO;
with Matrices;
procedure Multiple_Regression is
Line 318:
Output_Matrix (Float_Matrices.To_Matrix (Coefficients));
end;
end Multiple_Regression;</langsyntaxhighlight>
 
{{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|BBC BASIC}}==
{{works with|BBC BASIC for Windows}}
<langsyntaxhighlight lang="bbcbasic"> *FLOAT 64
INSTALL @lib$+"ARRAYLIB"
Line 389 ⟶ 605:
t() = t().m()
c() = y().t()
ENDPROC</langsyntaxhighlight>
{{out}}
<pre>
Line 396 ⟶ 612:
 
=={{header|C}}==
Using GNU gsl and c99, with the WP data<langsyntaxhighlight Clang="c">#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
Line 439 ⟶ 655:
gsl_multifit_linear_free(wspc);
 
}</langsyntaxhighlight>
 
=={{header|C++}}==
{{trans|Java}}
<syntaxhighlight lang="cpp">#include <array>
#include <iostream>
 
void require(bool condition, const std::string &message) {
if (condition) {
return;
}
throw std::runtime_error(message);
}
 
template<typename T, size_t N>
std::ostream &operator<<(std::ostream &os, const std::array<T, N> &a) {
auto it = a.cbegin();
auto end = a.cend();
 
os << '[';
if (it != end) {
os << *it;
it = std::next(it);
}
while (it != end) {
os << ", " << *it;
it = std::next(it);
}
return os << ']';
}
 
template <size_t RC, size_t CC>
class Matrix {
std::array<std::array<double, CC>, RC> data;
 
public:
Matrix() : data{} {
// empty
}
 
Matrix(std::initializer_list<std::initializer_list<double>> values) {
size_t rp = 0;
for (auto row : values) {
size_t cp = 0;
for (auto col : row) {
data[rp][cp] = col;
cp++;
}
rp++;
}
}
 
double get(size_t row, size_t col) const {
return data[row][col];
}
 
void set(size_t row, size_t col, double value) {
data[row][col] = value;
}
 
std::array<double, CC> get(size_t row) {
return data[row];
}
 
void set(size_t row, const std::array<double, CC> &values) {
std::copy(values.begin(), values.end(), data[row].begin());
}
 
template <size_t D>
Matrix<RC, D> operator*(const Matrix<CC, D> &rhs) const {
Matrix<RC, D> result;
for (size_t i = 0; i < RC; i++) {
for (size_t j = 0; j < D; j++) {
for (size_t k = 0; k < CC; k++) {
double prod = get(i, k) * rhs.get(k, j);
result.set(i, j, result.get(i, j) + prod);
}
}
}
return result;
}
 
Matrix<CC, RC> transpose() const {
Matrix<CC, RC> trans;
for (size_t i = 0; i < RC; i++) {
for (size_t j = 0; j < CC; j++) {
trans.set(j, i, data[i][j]);
}
}
return trans;
}
 
void toReducedRowEchelonForm() {
size_t lead = 0;
for (size_t r = 0; r < RC; r++) {
if (CC <= lead) {
return;
}
auto i = r;
 
while (get(i, lead) == 0.0) {
i++;
if (RC == i) {
i = r;
lead++;
if (CC == lead) {
return;
}
}
}
 
auto temp = get(i);
set(i, get(r));
set(r, temp);
 
if (get(r, lead) != 0.0) {
auto div = get(r, lead);
for (size_t j = 0; j < CC; j++) {
set(r, j, get(r, j) / div);
}
}
 
for (size_t k = 0; k < RC; k++) {
if (k != r) {
auto mult = get(k, lead);
for (size_t j = 0; j < CC; j++) {
auto prod = get(r, j) * mult;
set(k, j, get(k, j) - prod);
}
}
}
 
lead++;
}
}
 
Matrix<RC, RC> inverse() {
require(RC == CC, "Not a square matrix");
 
Matrix<RC, 2 * RC> aug;
for (size_t i = 0; i < RC; i++) {
for (size_t j = 0; j < RC; j++) {
aug.set(i, j, get(i, j));
}
// augment identify matrix to right
aug.set(i, i + RC, 1.0);
}
 
aug.toReducedRowEchelonForm();
 
// remove identity matrix to left
Matrix<RC, RC> inv;
for (size_t i = 0; i < RC; i++) {
for (size_t j = RC; j < 2 * RC; j++) {
inv.set(i, j - RC, aug.get(i, j));
}
}
return inv;
}
 
template <size_t RC, size_t CC>
friend std::ostream &operator<<(std::ostream &, const Matrix<RC, CC> &);
};
 
template <size_t RC, size_t CC>
std::ostream &operator<<(std::ostream &os, const Matrix<RC, CC> &m) {
for (size_t i = 0; i < RC; i++) {
os << '[';
for (size_t j = 0; j < CC; j++) {
if (j > 0) {
os << ", ";
}
os << m.get(i, j);
}
os << "]\n";
}
 
return os;
}
 
template <size_t RC, size_t CC>
std::array<double, RC> multiple_regression(const std::array<double, CC> &y, const Matrix<RC, CC> &x) {
Matrix<1, CC> tm;
tm.set(0, y);
 
auto cy = tm.transpose();
auto cx = x.transpose();
return ((x * cx).inverse() * x * cy).transpose().get(0);
}
 
void case1() {
std::array<double, 5> y{ 1.0, 2.0, 3.0, 4.0, 5.0 };
Matrix<1, 5> x{ {2.0, 1.0, 3.0, 4.0, 5.0} };
auto v = multiple_regression(y, x);
std::cout << v << '\n';
}
 
void case2() {
std::array<double, 3> y{ 3.0, 4.0, 5.0 };
Matrix<2, 3> x{
{1.0, 2.0, 1.0},
{1.0, 1.0, 2.0}
};
auto v = multiple_regression(y, x);
std::cout << v << '\n';
}
 
void case3() {
std::array<double, 15> 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 };
std::array<double, 15> a{ 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 };
 
Matrix<3, 15> x;
for (size_t i = 0; i < 15; i++) {
x.set(0, i, 1.0);
}
x.set(1, a);
for (size_t i = 0; i < 15; i++) {
x.set(2, i, a[i] * a[i]);
}
 
auto v = multiple_regression(y, x);
std::cout << v << '\n';
}
 
int main() {
case1();
case2();
case3();
 
return 0;
}</syntaxhighlight>
{{out}}
<pre>[0.981818]
[1, 2]
[128.813, -143.162, 61.9603]</pre>
 
=={{header|C sharp|C#}}==
{{libheader|Math.Net}}
<langsyntaxhighlight lang="csharp">using System;
using MathNet.Numerics.LinearRegression;
using MathNet.Numerics.LinearAlgebra;
Line 460 ⟶ 910:
Console.WriteLine(β);
}
}</langsyntaxhighlight>
 
{{out}}
Line 467 ⟶ 917:
-143.162
61.9603</pre>
 
 
=={{header|Common Lisp}}==
Line 473 ⟶ 922:
Uses the routine (chol A) from [[Cholesky decomposition]], (mmul A B) from [[Matrix multiplication]], (mtp A) from [[Matrix transposition]].
 
<langsyntaxhighlight lang="lisp">
;; Solve a linear system AX=B where A is symmetric and positive definite, so it can be Cholesky decomposed.
(defun linsys (A B)
Line 507 ⟶ 956:
(linsys (mmul (mtp A) A)
(mmul (mtp A) b)))
</syntaxhighlight>
</lang>
 
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.
 
<langsyntaxhighlight lang="lisp">(let ((x (make-array '(1 11) :initial-contents '((0 1 2 3 4 5 6 7 8 9 10))))
(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))</langsyntaxhighlight>
 
=={{header|Emacs LispD}}==
{{trans|Java}}
<syntaxhighlight lang="d">import std.algorithm;
import std.array;
import std.exception;
import std.range;
import std.stdio;
 
public class Matrix {
Multiple regression analysis by Emacs Lisp and built-in Emacs Calc.
private double[][] data;
private size_t rowCount;
private size_t colCount;
 
public this(size_t size)
<lang emacs-lisp>
(setq X1 '[0 1 2in(size 3> 40, 5"Must 6have 7at 8least 9one 10]element")
{
(setq X2 '[0 1 1 3 3 7 6 7 3 9 8])
(setq Y '[1 6 17 34 57 86 121 162 209 262this(size, 321]size);
}
(calc-eval
(format "fit(a*X1+b*X2+c,[X1,X2],[a,b,c],[%s %s %s])" X1 X2 Y))
</lang>
 
public this(size_t rows, size_t cols)
in(rows > 0, "Must have at least one row")
in(cols > 0, "Must have at least one column")
{
rowCount = rows;
colCount = cols;
 
data = uninitializedArray!(double[][])(rows, cols);
foreach (ref row; data) {
row[] = 0.0;
}
}
 
public this(const double[][] source) {
enforce(source.length > 0, "Must have at least one row");
rowCount = source.length;
 
enforce(source[0].length > 0, "Must have at least one column");
colCount = source[0].length;
 
data = uninitializedArray!(double[][])(rowCount, colCount);
foreach (i; 0 .. rowCount) {
enforce(source[i].length == colCount, "All rows must have equal columns");
data[i] = source[i].dup;
}
}
 
public auto opIndex(size_t r, size_t c) const {
return data[r][c];
}
 
public auto opIndex(size_t r) const {
return data[r];
}
 
public auto opBinary(string op)(const Matrix rhs) const {
static if (op == "*") {
auto rc1 = rowCount;
auto cc1 = colCount;
auto rc2 = rhs.rowCount;
auto cc2 = rhs.colCount;
enforce(cc1 == rc2, "Cannot multiply if the first columns does not equal the second rows");
auto result = new Matrix(rc1, cc2);
foreach (i; 0 .. rc1) {
foreach (j; 0 .. cc2) {
foreach (k; 0 .. rc2) {
result[i, j] += this[i, k] * rhs[k, j];
}
}
}
return result;
} else {
assert(false, "Not implemented");
}
}
 
public void opIndexAssign(double value, size_t r, size_t c) {
data[r][c] = value;
}
 
public void opIndexAssign(const double[] value, size_t r) {
enforce(colCount == value.length, "Slice size must match column size");
data[r] = value.dup;
}
 
public void opIndexOpAssign(string op)(double value, size_t r, size_t c) {
mixin("data[r][c] " ~ op ~ "= value;");
}
 
public auto transpose() const {
auto rc = rowCount;
auto cc = colCount;
auto t = new Matrix(cc, rc);
foreach (i; 0 .. cc) {
foreach (j; 0 .. rc) {
t[i, j] = this[j, i];
}
}
return t;
}
 
public void toReducedRowEchelonForm() {
auto lead = 0;
auto rc = rowCount;
auto cc = colCount;
foreach (r; 0 .. rc) {
if (cc <= lead) {
return;
}
auto i = r;
 
while (this[i, lead] == 0.0) {
i++;
if (rc == i) {
i = r;
lead++;
if (cc == lead) {
return;
}
}
}
 
auto temp = this[i];
this[i] = this[r];
this[r] = temp;
 
if (this[r, lead] != 0.0) {
auto div = this[r, lead];
foreach (j; 0 .. cc) {
this[r, j] = this[r, j] / div;
}
}
 
foreach (k; 0 .. rc) {
if (k != r) {
auto mult = this[k, lead];
foreach (j; 0 .. cc) {
this[k, j] -= this[r, j] * mult;
}
}
}
 
lead++;
}
}
 
public auto inverse() const {
enforce(rowCount == colCount, "Not a square matrix");
auto len = rowCount;
auto aug = new Matrix(len, 2 * len);
foreach (i; 0 .. len) {
foreach (j; 0 .. len) {
aug[i, j] = this[i, j];
}
// augment identity matrix to right
aug[i, i + len] = 1.0;
}
aug.toReducedRowEchelonForm;
auto inv = new Matrix(len);
// remove identify matrix to left
foreach (i; 0 .. len) {
foreach (j; len .. 2 * len) {
inv[i, j - len] = aug[i, j];
}
}
return inv;
}
 
void toString(scope void delegate(const(char)[]) sink) const {
import std.format;
auto fmt = FormatSpec!char("%s");
 
put(sink, "[");
foreach (i; 0 .. rowCount) {
if (i > 0) {
put(sink, " [");
} else {
put(sink, "[");
}
 
formatValue(sink, this[i, 0], fmt);
foreach (j; 1 .. colCount) {
put(sink, ", ");
formatValue(sink, this[i, j], fmt);
}
 
if (i + 1 < rowCount) {
put(sink, "]\n");
} else {
put(sink, "]");
}
}
put(sink, "]");
}
}
 
auto multipleRegression(double[] y, Matrix x) {
auto tm = new Matrix([y]);
auto cy = tm.transpose;
auto cx = x.transpose;
return ((x * cx).inverse * x * cy).transpose[0].dup;
}
 
void main() {
auto y = [1.0, 2.0, 3.0, 4.0, 5.0];
auto x = new Matrix([[2.0, 1.0, 3.0, 4.0, 5.0]]);
auto v = multipleRegression(y, x);
v.writeln;
 
y = [3.0, 4.0, 5.0];
x = new Matrix([
[1.0, 2.0, 1.0],
[1.0, 1.0, 2.0]
]);
v = multipleRegression(y, x);
v.writeln;
 
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];
auto a = [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];
x = new Matrix([
repeat(1.0, a.length).array,
a,
a.map!"a * a".array
]);
v = multipleRegression(y, x);
v.writeln;
}</syntaxhighlight>
{{out}}
<pre>[0.981818]
[1, 2]
"35.2014388489 X1 - 3.95683453237 X2 - 42.7410071942"
[128.813, -143.162, 61.9603]</pre>
</pre>
 
=={{header|Emacs Lisp}}==
 
{{libheader|calc}}
 
<syntaxhighlight lang="lisp">(let ((x1 '(0 1 2 3 4 5 6 7 8 9 10))
(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))))</syntaxhighlight>
 
{{out}}
 
"35.2014388489*X1 - 3.95683453237*X2 - 42.7410071942"
 
=={{header|ERRE}}==
<langsyntaxhighlight ERRElang="erre">PROGRAM MULTIPLE_REGRESSION
 
!$DOUBLE
Line 618 ⟶ 1,293:
END FOR
 
END PROGRAM</langsyntaxhighlight>
{{out}}
<pre>LINEAR SYSTEM COEFFICENTS
Line 638 ⟶ 1,313:
{{libheader|SLATEC}} [http://netlib.org/slatec/ Available at the Netlib]
 
<langsyntaxhighlight Fortranlang="fortran">*-----------------------------------------------------------------------
* MR - multiple regression using the SLATEC library routine DHFTI
*
Line 723 ⟶ 1,398:
STOP 'program complete'
END
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 729 ⟶ 1,404:
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===
<langsyntaxhighlight lang="go">package main
 
import (
Line 765 ⟶ 1,517:
x, y := givens()
fmt.Printf("%.4f\n", mat64.Formatted(mat64.QR(x).Solve(y)))
}</langsyntaxhighlight>
{{out}}
<pre>
Line 773 ⟶ 1,525:
</pre>
===Library go.matrix===
<langsyntaxhighlight lang="go">package main
 
import (
Line 818 ⟶ 1,570:
}
fmt.Println(c)
}</langsyntaxhighlight>
{{out}}
<pre>
Line 826 ⟶ 1,578:
=={{header|Haskell}}==
Using package [http://hackage.haskell.org/package/hmatrix hmatrix] from HackageDB
<langsyntaxhighlight lang="haskell">import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.LAPACK
 
Line 837 ⟶ 1,589:
v :: Matrix Double
v = (3><1)
[1.745005,-4.448092,-4.160842]</langsyntaxhighlight>
Using lapack::dgels
<langsyntaxhighlight lang="haskell">*Main> linearSolveLSR m v
(3><1)
[ 0.9335611922087276
, 1.101323491272865
, 1.6117769115824 ]</langsyntaxhighlight>
Or
<langsyntaxhighlight lang="haskell">*Main> inv m `multiply` v
(3><1)
[ 0.9335611922087278
, 1.101323491272865
, 1.6117769115824006 ]</langsyntaxhighlight>
 
=={{header|Hy}}==
<langsyntaxhighlight lang="lisp">(import
[numpy [ones column-stack]]
[numpy.random [randn]]
Line 864 ⟶ 1,616:
(print (first (lstsq
(column-stack (, (ones n) x1 x2 (* x1 x2)))
y)))</langsyntaxhighlight>
 
=={{header|J}}==
 
<langsyntaxhighlight lang="j"> NB. Wikipedia data
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</langsyntaxhighlight>
 
Breaking it down:
<langsyntaxhighlight lang="j"> X=: x ^/ i.3 NB. form Design matrix
X=: (x^0) ,. (x^1) ,. (x^2) NB. equivalent of previous line
4{.X NB. show first 4 rows of X
Line 887 ⟶ 1,639:
NB. y %. X does matrix division and gives the regression coefficients
y %. X
128.813 _143.162 61.9603</langsyntaxhighlight>
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:
<langsyntaxhighlight lang="j"> mp=: +/ .* NB. matrix product
NB. %.X is matrix inverse of X
NB. |:X is transpose of X
Line 901 ⟶ 1,653:
X (%.@:xpy@[ mp xpy) y
128.814 _143.163 61.9606
</syntaxhighlight>
</lang>
 
LAPACK routines are also available via the Addon <tt>math/lapack</tt>.
<langsyntaxhighlight lang="j"> load 'math/lapack'
load 'math/lapack/gels'
gels_jlapack_ X;y
128.813 _143.162 61.9603</langsyntaxhighlight>
 
=={{header|JuliaJava}}==
{{trans|MATLABKotlin}}
<syntaxhighlight lang="java">import java.util.Arrays;
import java.util.Objects;
 
public class MultipleRegression {
As in Matlab, the backslash or slash operator (depending on the matrix ordering) can be used for solving this problem, for example:
public static void require(boolean condition, String message) {
if (condition) {
return;
}
throw new IllegalArgumentException(message);
}
 
public static class Matrix {
<lang julia>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]
private final double[][] data;
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]
private final int rowCount;
X = [x.^0 x.^1 x.^2];
private final int colCount;
b = X \ y</lang>
 
public Matrix(int rows, int cols) {
require(rows > 0, "Need at least one row");
this.rowCount = rows;
 
require(cols > 0, "Need at least one column");
this.colCount = cols;
 
this.data = new double[rows][cols];
for (double[] row : this.data) {
Arrays.fill(row, 0.0);
}
}
 
public Matrix(double[][] source) {
require(source.length > 0, "Need at least one row");
this.rowCount = source.length;
 
require(source[0].length > 0, "Need at least one column");
this.colCount = source[0].length;
 
this.data = new double[this.rowCount][this.colCount];
for (int i = 0; i < this.rowCount; i++) {
set(i, source[i]);
}
}
 
public double[] get(int row) {
Objects.checkIndex(row, this.rowCount);
return this.data[row];
}
 
public void set(int row, double[] data) {
Objects.checkIndex(row, this.rowCount);
require(data.length == this.colCount, "The column in the row must match the number of columns in the matrix");
System.arraycopy(data, 0, this.data[row], 0, this.colCount);
}
 
public double get(int row, int col) {
Objects.checkIndex(row, this.rowCount);
Objects.checkIndex(col, this.colCount);
return this.data[row][col];
}
 
public void set(int row, int col, double value) {
Objects.checkIndex(row, this.rowCount);
Objects.checkIndex(col, this.colCount);
this.data[row][col] = value;
}
 
@SuppressWarnings("UnnecessaryLocalVariable")
public Matrix times(Matrix that) {
var rc1 = this.rowCount;
var cc1 = this.colCount;
var rc2 = that.rowCount;
var cc2 = that.colCount;
require(cc1 == rc2, "Cannot multiply if the first columns does not equal the second rows");
var result = new Matrix(rc1, cc2);
for (int i = 0; i < rc1; i++) {
for (int j = 0; j < cc2; j++) {
for (int k = 0; k < rc2; k++) {
var prod = get(i, k) * that.get(k, j);
result.set(i, j, result.get(i, j) + prod);
}
}
}
return result;
}
 
public Matrix transpose() {
var rc = this.rowCount;
var cc = this.colCount;
var trans = new Matrix(cc, rc);
for (int i = 0; i < cc; i++) {
for (int j = 0; j < rc; j++) {
trans.set(i, j, get(j, i));
}
}
return trans;
}
 
public void toReducedRowEchelonForm() {
int lead = 0;
var rc = this.rowCount;
var cc = this.colCount;
for (int r = 0; r < rc; r++) {
if (cc <= lead) {
return;
}
var i = r;
 
while (get(i, lead) == 0.0) {
i++;
if (rc == i) {
i = r;
lead++;
if (cc == lead) {
return;
}
}
}
 
var temp = get(i);
set(i, get(r));
set(r, temp);
 
if (get(r, lead) != 0.0) {
var div = get(r, lead);
for (int j = 0; j < cc; j++) {
set(r, j, get(r, j) / div);
}
}
 
for (int k = 0; k < rc; k++) {
if (k != r) {
var mult = get(k, lead);
for (int j = 0; j < cc; j++) {
var prod = get(r, j) * mult;
set(k, j, get(k, j) - prod);
}
}
}
 
lead++;
}
}
 
public Matrix inverse() {
require(this.rowCount == this.colCount, "Not a square matrix");
var len = this.rowCount;
var aug = new Matrix(len, 2 * len);
for (int i = 0; i < len; i++) {
for (int j = 0; j < len; j++) {
aug.set(i, j, get(i, j));
}
// augment identity matrix to right
aug.set(i, i + len, 1.0);
}
aug.toReducedRowEchelonForm();
var inv = new Matrix(len, len);
// remove identity matrix to left
for (int i = 0; i < len; i++) {
for (int j = len; j < 2 * len; j++) {
inv.set(i, j - len, aug.get(i, j));
}
}
return inv;
}
}
 
public static double[] multipleRegression(double[] y, Matrix x) {
var tm = new Matrix(new double[][]{y});
var cy = tm.transpose();
var cx = x.transpose();
return x.times(cx).inverse().times(x).times(cy).transpose().get(0);
}
 
public static void printVector(double[] v) {
System.out.println(Arrays.toString(v));
System.out.println();
}
 
public static double[] repeat(int size, double value) {
var a = new double[size];
Arrays.fill(a, value);
return a;
}
 
public static void main(String[] args) {
double[] y = new double[]{1.0, 2.0, 3.0, 4.0, 5.0};
var x = new Matrix(new double[][]{{2.0, 1.0, 3.0, 4.0, 5.0}});
var v = multipleRegression(y, x);
printVector(v);
 
y = new double[]{3.0, 4.0, 5.0};
x = new Matrix(new double[][]{
{1.0, 2.0, 1.0},
{1.0, 1.0, 2.0}
});
v = multipleRegression(y, x);
printVector(v);
 
y = new double[]{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};
var a = new double[]{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};
x = new Matrix(new double[][]{
repeat(a.length, 1.0),
a,
Arrays.stream(a).map(it -> it * it).toArray()
});
 
v = multipleRegression(y, x);
printVector(v);
}
}</syntaxhighlight>
{{out}}
<pre>[0.9818181818181818]
<pre>
 
3-element Array{Float64,1}:
[0.9999999999999996, 2.000000000000001]
128.813
 
-143.162
[128.8128035798277, -143.1620228653037, 61.960325442985436]</pre>
61.9603
</pre>
 
=={{header|JavaScript}}==
Line 933 ⟶ 1,886:
Extends the Matrix class from [[Matrix Transpose#JavaScript]], [[Matrix multiplication#JavaScript]], [[Reduced row echelon form#JavaScript]].
Uses the IdentityMatrix from [[Matrix exponentiation operator#JavaScript]]
<langsyntaxhighlight lang="javascript">// modifies the matrix "in place"
Matrix.prototype.inverse = function() {
if (this.height != this.width) {
Line 979 ⟶ 1,932:
)
);
print(y.regression_coefficients(x));</langsyntaxhighlight>
{{out}}
<pre>0.9818181818181818
Line 987 ⟶ 1,940:
61.960325442985436</pre>
 
=={{header|Mathematicajq}}==
{{trans|Wren}}
<lang Mathematica>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};
{{works with|jq}}
'''Works with gojq, the Go implementation of jq'''
 
See e.g. https://rosettacode.org/wiki/Gauss-Jordan_matrix_inversion#jq
for a suitable definition of `inverse` as used here.
 
'''Preliminaries'''
<syntaxhighlight lang="jq">def dot_product(a; b):
reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );
 
# A and B should both be numeric matrices, A being m by n, and B being n by p.
def multiply(A; B):
(B[0]|length) as $p
| (B|transpose) as $BT
| reduce range(0; A|length) as $i
([];
reduce range(0; $p) as $j
(.;
.[$i][$j] = dot_product( A[$i]; $BT[$j] ) ));</syntaxhighlight>
 
'''Multiple Regression'''
<syntaxhighlight lang="jq">def multipleRegression(y; x):
(y|transpose) as $cy
| (x|transpose) as $cx
| multiply( multiply( multiply(x;$cx)|inverse; x); $cy)
| transpose[0];
def ys: [
[ [1, 2, 3, 4, 5] ],
[ [3, 4, 5] ],
[ [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] ]
];
def a:
[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];
 
def xs:[
[ [2, 1, 3, 4, 5] ],
[ [1, 2, 1], [1, 1, 2] ],
[ [range(0;a|length) | 1], a, (a|map(.*.))]
];
 
range(0; ys|length) as $i
| multipleRegression(ys[$i]; xs[$i])</syntaxhighlight>
{{out}}
<pre>
[0.9818181818181818]
[0.9999999999999996,2.000000000000001]
[128.8128035798277,-143.1620228653037,61.960325442985436]
</pre>
 
 
=={{header|Julia}}==
{{trans|MATLAB}}
 
As in Matlab, the backslash or slash operator (depending on the matrix ordering) can be used for solving this problem, for example:
 
<syntaxhighlight lang="julia">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]
X = [x.^0 x.^1 x.^2];
b = X \ y</syntaxhighlight>
{{out}}
<pre>
3-element Array{Float64,1}:
128.813
-143.162
61.9603
</pre>
 
=={{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.
<syntaxhighlight lang="scala">// Version 1.2.31
 
typealias Vector = DoubleArray
typealias Matrix = Array<Vector>
operator fun Matrix.times(other: Matrix): Matrix {
val rows1 = this.size
val cols1 = this[0].size
val rows2 = other.size
val cols2 = other[0].size
require(cols1 == rows2)
val result = Matrix(rows1) { Vector(cols2) }
for (i in 0 until rows1) {
for (j in 0 until cols2) {
for (k in 0 until rows2) {
result[i][j] += this[i][k] * other[k][j]
}
}
}
return result
}
 
fun Matrix.transpose(): Matrix {
val rows = this.size
val cols = this[0].size
val trans = Matrix(cols) { Vector(rows) }
for (i in 0 until cols) {
for (j in 0 until rows) trans[i][j] = this[j][i]
}
return trans
}
 
fun Matrix.inverse(): Matrix {
val len = this.size
require(this.all { it.size == len }) { "Not a square matrix" }
val aug = Array(len) { DoubleArray(2 * len) }
for (i in 0 until len) {
for (j in 0 until len) aug[i][j] = this[i][j]
// augment by identity matrix to right
aug[i][i + len] = 1.0
}
aug.toReducedRowEchelonForm()
val inv = Array(len) { DoubleArray(len) }
// remove identity matrix to left
for (i in 0 until len) {
for (j in len until 2 * len) inv[i][j - len] = aug[i][j]
}
return inv
}
 
fun Matrix.toReducedRowEchelonForm() {
var lead = 0
val rowCount = this.size
val colCount = this[0].size
for (r in 0 until rowCount) {
if (colCount <= lead) return
var i = r
 
while (this[i][lead] == 0.0) {
i++
if (rowCount == i) {
i = r
lead++
if (colCount == lead) return
}
}
 
val temp = this[i]
this[i] = this[r]
this[r] = temp
 
if (this[r][lead] != 0.0) {
val div = this[r][lead]
for (j in 0 until colCount) this[r][j] /= div
}
for (k in 0 until rowCount) {
if (k != r) {
val mult = this[k][lead]
for (j in 0 until colCount) this[k][j] -= this[r][j] * mult
}
}
 
lead++
}
}
 
fun printVector(v: Vector) {
println(v.asList())
println()
}
 
fun multipleRegression(y: Vector, x: Matrix): Vector {
val cy = (arrayOf(y)).transpose() // convert 'y' to column vector
val cx = x.transpose() // convert 'x' to column vector array
return ((x * cx).inverse() * x * cy).transpose()[0]
}
 
fun main(args: Array<String>) {
var y = doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0)
var x = arrayOf(doubleArrayOf(2.0, 1.0, 3.0, 4.0, 5.0))
var v = multipleRegression(y, x)
printVector(v)
 
y = doubleArrayOf(3.0, 4.0, 5.0)
x = arrayOf(
doubleArrayOf(1.0, 2.0, 1.0),
doubleArrayOf(1.0, 1.0, 2.0)
)
v = multipleRegression(y, x)
printVector(v)
 
y = doubleArrayOf(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)
 
val a = doubleArrayOf(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)
x = arrayOf(DoubleArray(a.size) { 1.0 }, a, a.map { it * it }.toDoubleArray())
v = multipleRegression(y, x)
printVector(v)
}</syntaxhighlight>
 
{{out}}
<pre>
[0.9818181818181818]
 
[0.9999999999999996, 2.000000000000001]
 
[128.8128035798277, -143.1620228653037, 61.960325442985436]
</pre>
 
=={{header|Maple}}==
 
First build a random dataset.
 
<syntaxhighlight lang="maple">n:=200:
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):</syntaxhighlight>
 
Now the linear regression, with either the LinearAlgebra package, or the Statistics package.
 
<syntaxhighlight lang="maple">LinearAlgebra[LeastSquares](X,Y)^+;
# [4.33701132468683, -2.78654498997457, -1.41840666085642, 2.92065133466547, 1.76076442997642]
 
Statistics[LinearFit]([x1,x2,x3,x4,c],X,Y,[x1,x2,x3,x4,c],summarize=true)
# Summary:
# ----------------
# Model: 4.3370113*x1-2.7865450*x2-1.4184067*x3+2.9206513*x4+1.7607644*c
# ----------------
# Coefficients:
# Estimate Std. Error t-value P(>|t|)
# Parameter 1 4.3370 0.0691 62.7409 0.0000
# Parameter 2 -2.7865 0.0661 -42.1637 0.0000
# Parameter 3 -1.4184 0.0699 -20.2937 0.0000
# Parameter 4 2.9207 0.0687 42.5380 0.0000
# Parameter 5 1.7608 0.0701 25.1210 0.0000
# ----------------
# R-squared: 0.9767, Adjusted R-squared: 0.9761
# 4.33701132468683 x1 - 2.78654498997457 x2 - 1.41840666085642 x3
# + 2.92065133466547 x4 + 1.76076442997642 c</syntaxhighlight>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">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};
X = {x^0, x^1, x^2};
LeastSquares[Transpose@X, y]</syntaxhighlight>
b = y.PseudoInverse[X]
{{out}}
 
-<pre>{128.813, -143.162, 61.9603}</langpre>
 
=={{header|MATLAB}}==
Line 999 ⟶ 2,187:
The slash and backslash operator can be used for solving this problem. Here some random data is generated.
 
<langsyntaxhighlight Matlablang="matlab"> n=100; k=10;
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 </langsyntaxhighlight>
 
In its transposed form yt = Xt * bt, the backslash operator can be used.
 
<langsyntaxhighlight Matlablang="matlab"> yt = y'; Xt = X';
bt = Xt \ yt
bt =
Line 1,019 ⟶ 2,207:
-0.0561589
-0.1752146
-0.2577663</langsyntaxhighlight>
 
Here is the example for estimating the polynomial fit
 
<langsyntaxhighlight Matlablang="matlab"> 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]
X = [x.^0;x.^1;x.^2];
b = y/X
 
128.813 -143.162 61.960</langsyntaxhighlight>
 
Instead of "/", the slash operator, one can also write :
<langsyntaxhighlight Matlablang="matlab"> b = y * X' * inv(X * X') </langsyntaxhighlight>
or
<langsyntaxhighlight Matlablang="matlab"> b = y * pinv(X) </langsyntaxhighlight>
 
=={{header|Nim}}==
{{libheader|arraymancer}}
<syntaxhighlight lang="nim"># Using Wikipedia data sample.
 
import math
import arraymancer, sequtils
 
var
 
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].toTensor()
 
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].toTensor()
 
# Create Vandermonde matrix.
var a = stack(height.ones_like, height, height *. height, axis = 1)
 
echo toSeq(least_squares_solver(a, weight).solution.items)</syntaxhighlight>
 
{{out}}
<pre>@[128.8128035784397, -143.1620228647656, 61.96032544247468]</pre>
 
=={{header|PARI/GP}}==
<langsyntaxhighlight lang="parigp">pseudoinv(M)=my(sz=matsize(M),T=conj(M))~;if(sz[1]<sz[2],T/(M*T),(T*M)^-1*T)
addhelp(pseudoinv, "pseudoinv(M): Moore pseudoinverse of the matrix M.");
 
y*pseudoinv(X)</langsyntaxhighlight>
 
=={{header|Perl}}==
<syntaxhighlight lang="perl">use strict;
use warnings;
use Statistics::Regression;
 
my @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);
my @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);
 
my @model = ('const', 'X', 'X**2');
my $reg = Statistics::Regression->new( '', [@model] );
$reg->include( $y[$_], [ 1.0, $x[$_], $x[$_]**2 ]) for 0..@y-1;
my @coeff = $reg->theta();
 
printf "%-6s %8.3f\n", $model[$_], $coeff[$_] for 0..@model-1;</syntaxhighlight>
{{out}}
<pre>const 128.813
X -143.162
X**2 61.960</pre>
 
=={{header|Phix}}==
{{trans|ERRE}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<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>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1.47</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.50</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.52</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.55</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.57</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">1.60</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.63</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.65</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.68</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.70</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">1.73</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.75</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.78</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.80</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1.83</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">52.21</span><span style="color: #0000FF;">,</span><span style="color: #000000;">53.12</span><span style="color: #0000FF;">,</span><span style="color: #000000;">54.48</span><span style="color: #0000FF;">,</span><span style="color: #000000;">55.84</span><span style="color: #0000FF;">,</span><span style="color: #000000;">57.20</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">58.57</span><span style="color: #0000FF;">,</span><span style="color: #000000;">59.93</span><span style="color: #0000FF;">,</span><span style="color: #000000;">61.29</span><span style="color: #0000FF;">,</span><span style="color: #000000;">63.11</span><span style="color: #0000FF;">,</span><span style="color: #000000;">64.47</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">66.28</span><span style="color: #0000FF;">,</span><span style="color: #000000;">68.10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">69.92</span><span style="color: #0000FF;">,</span><span style="color: #000000;">72.19</span><span style="color: #0000FF;">,</span><span style="color: #000000;">74.46</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">N</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">N</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</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;">M</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">M</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">N</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">k</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">M</span> <span style="color: #008080;">then</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]*</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000080;font-style:italic;">-- build linear system</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">row</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">M</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">col</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">M</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">row</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">row</span><span style="color: #0000FF;">+</span><span style="color: #000000;">col</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">row</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: #0000FF;">=</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">[</span><span style="color: #000000;">row</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<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;">"Linear system coefficents:\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,{</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%7.1f"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_FltFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%7.1f"</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">M</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">j</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span> <span style="color: #000000;">i</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">i</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: #008080;">then</span>
<span style="color: #0000FF;">?</span><span style="color: #008000;">"SINGULAR MATRIX !"</span>
<span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">M</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">Y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span><span style="color: #000000;">Y</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">M</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">k</span><span style="color: #0000FF;"><></span><span style="color: #000000;">j</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">Y</span><span style="color: #0000FF;">=-</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">M</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">Y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<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>
<!--</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.96032544}
</pre>
 
=={{header|PicoLisp}}==
<langsyntaxhighlight PicoLisplang="picolisp">(scl 20)
 
# Matrix transposition
Line 1,084 ⟶ 2,383:
(car X) ) ) ) )
(T (> (inc 'Lead) Cols)) ) )
Mat )</langsyntaxhighlight>
{{trans|JavaScript}}
<langsyntaxhighlight PicoLisplang="picolisp">(de matInverse (Mat)
(let N (length Mat)
(unless (= N (length (car Mat)))
Line 1,104 ⟶ 2,403:
X (columnVector (2.0 1.0 3.0 4.0 5.0)) )
 
(round (caar (regressionCoefficients Y X)) 17)</langsyntaxhighlight>
{{out}}
<pre>-> "0.98181818181818182"</pre>
 
 
=={{header|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.
 
Let's create four vectors containing our input data:
 
<math>\begin{align}
\mathbf{w} & = w^k\mathbf{e}_k\\
\mathbf{h_0} & = (h^k)^0\mathbf{e}_k\\
\mathbf{h_1} & = (h^k)^1\mathbf{e}_k\\
\mathbf{h_2} & = (h^k)^2\mathbf{e}_k
\end{align}</math>
 
Then what we're looking for are three scalars <math>\alpha</math>, <math>\beta</math> and <math>\gamma</math> such that:
 
<math>\alpha\mathbf{h0} + \beta\mathbf{h1} + \gamma\mathbf{h2} = \mathbf{w}</math>
 
To get for instance <math>\alpha</math> we can first make the <math>\beta</math> and <math>\gamma</math> terms disappear:
 
<math>\alpha\mathbf{h0}\wedge\mathbf{h1}\wedge\mathbf{h2} = \mathbf{w}\wedge\mathbf{h1}\wedge\mathbf{h2}</math>
 
Noting <math>I = \mathbf{h0}\wedge\mathbf{h1}\wedge\mathbf{h2}</math>, we then get:
 
<math>\alpha = (\mathbf{w}\wedge\mathbf{h1}\wedge\mathbf{h2})\cdot\tilde{I}/I\cdot\tilde{I}</math>
 
'''Note:''' a number of the formulae above are invisible to the majority of browsers, including Chrome, IE/Edge, Safari and Opera. They may (subject to the installation of necessary fronts) be visible to Firefox.
 
 
<lang perl6>use Clifford;
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>;
 
my $w = [+] @weight Z* @e;
 
my $h0 = [+] @e[^@weight];
my $h1 = [+] @height Z* @e;
my $h2 = [+] (@height X** 2) Z* @e;
 
my $I = $h0∧$h1∧$h2;
my $I2 = ($I·$I.reversion).Real;
 
say "α = ", ($w∧$h1∧$h2)·$I.reversion/$I2;
say "β = ", ($w∧$h2∧$h0)·$I.reversion/$I2;
say "γ = ", ($w∧$h0∧$h1)·$I.reversion/$I2;</lang>
{{out}}
<pre>α = 128.81280357844
β = -143.1620228648
γ = 61.960325442</pre>
 
This computation took over an hour with the april 2016 version of rakudo on MoarVM, running in a VirtualBox linux system guest hosted by a windows laptop with a i7 intel processor.
 
=={{header|Python}}==
{{libheader|NumPy}}
'''Method with matrix operations'''
<langsyntaxhighlight lang="python">import numpy as np
 
height = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,
Line 1,172 ⟶ 2,420:
y = np.mat(weight)
 
print(y * X.T * (X*X.T).I)</langsyntaxhighlight>
{{out}}
<pre>
Line 1,178 ⟶ 2,426:
</pre>
'''Using numpy lstsq function'''
<langsyntaxhighlight lang="python">import numpy as np
 
height = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,
Line 1,188 ⟶ 2,436:
y = weight
 
print(np.linalg.lstsq(X, y)[0])</langsyntaxhighlight>
{{out}}
<pre>
Line 1,196 ⟶ 2,444:
=={{header|R}}==
 
R provides the '''lm()''' function for linear regression.
 
<syntaxhighlight lang="rsplus">x <- c(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)
<lang R>## Wikipedia Data
x <- c(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 <- 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))</langsyntaxhighlight>
 
{{out}}
 
<pre>
<pre>Call:
lm(formula = y ~ x + I(x^2))
 
Coefficients:
(Intercept) x I(x^2)
128.81 -143.16 61.96 </pre>
</pre>
 
A simple implementation of multiple regression in native R
is useful to illustrate R's model description and linear algebra capabilities.
 
<langsyntaxhighlight Rlang="rsplus">simpleMultipleReg <- function(formula) {
 
## parse and evaluate the model formula
Line 1,223 ⟶ 2,469:
 
## create design matrix
X <- model.matrix(attr(mf, "terms"), mf)
 
## create dependent variable
Line 1,232 ⟶ 2,478:
}
 
simpleMultipleReg(y ~ x + I(x^2))</langsyntaxhighlight>
 
This produces the same coefficients as lm()
Line 1,246 ⟶ 2,492:
than the method above, is to solve the linear system directly
and use the crossprod function:
<langsyntaxhighlight Rlang="r">solve( crossprod(X), crossprod(X, Y))</langsyntaxhighlight>
 
A numerically more stable way is to use the QR decomposition of the design matrix:
 
<syntaxhighlight lang="rsplus">lm.impl <- function(formula) {
mf <- model.frame(formula)
X <- model.matrix(mf)
Y <- model.response(mf)
qr.coef(qr(X), Y)
}
 
 
lm(y ~ x + I(x^2))
 
# Call:
# lm(formula = y ~ x + I(x^2))
#
# Coefficients:
# (Intercept) x I(x^2)
# 128.81 -143.16 61.96
 
lm.impl(y ~ x + I(x^2))
 
# (Intercept) x I(x^2)
# 128.81280 -143.16202 61.96033</syntaxhighlight>
 
=={{header|Racket}}==
<langsyntaxhighlight lang="racket">
#lang racket
(require math)
Line 1,256 ⟶ 2,526:
(define (fit X y)
(matrix-solve (matrix* (T X) X) (matrix* (T X) y)))
</syntaxhighlight>
</lang>
Test:
<langsyntaxhighlight lang="racket">
(fit (matrix [[1 2]
[2 5]
Line 1,269 ⟶ 2,539:
{{out}}
(array #[#[9 1/3] #[-3 1/3]])
</syntaxhighlight>
</lang>
 
=={{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.
 
Let's create four vectors containing our input data:
 
<math>\begin{align}
\mathbf{w} & = w^k\mathbf{e}_k\\
\mathbf{h_0} & = (h^k)^0\mathbf{e}_k\\
\mathbf{h_1} & = (h^k)^1\mathbf{e}_k\\
\mathbf{h_2} & = (h^k)^2\mathbf{e}_k
\end{align}</math>
 
Then what we're looking for are three scalars <math>\alpha</math>, <math>\beta</math> and <math>\gamma</math> such that:
 
<math>\alpha\mathbf{h0} + \beta\mathbf{h1} + \gamma\mathbf{h2} = \mathbf{w}</math>
 
To get for instance <math>\alpha</math> we can first make the <math>\beta</math> and <math>\gamma</math> terms disappear:
 
<math>\alpha\mathbf{h0}\wedge\mathbf{h1}\wedge\mathbf{h2} = \mathbf{w}\wedge\mathbf{h1}\wedge\mathbf{h2}</math>
 
Noting <math>I = \mathbf{h0}\wedge\mathbf{h1}\wedge\mathbf{h2}</math>, we then get:
 
<math>\alpha = (\mathbf{w}\wedge\mathbf{h1}\wedge\mathbf{h2})\cdot\tilde{I}/I\cdot\tilde{I}</math>
 
'''Note:''' a number of the formulae above are invisible to the majority of browsers, including Chrome, IE/Edge, Safari and Opera. They may (subject to the installation of necessary fronts) be visible to Firefox.
 
 
<syntaxhighlight lang="raku" line>use Clifford;
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>;
 
my $w = [+] @weight Z* @e;
 
my $h0 = [+] @e[^@weight];
my $h1 = [+] @height Z* @e;
my $h2 = [+] (@height X** 2) Z* @e;
 
my $I = $h0∧$h1∧$h2;
my $I2 = ($I·$I.reversion).Real;
 
say "α = ", ($w∧$h1∧$h2)·$I.reversion/$I2;
say "β = ", ($w∧$h2∧$h0)·$I.reversion/$I2;
say "γ = ", ($w∧$h0∧$h1)·$I.reversion/$I2;</syntaxhighlight>
{{out}}
<pre>α = 128.81280357844
β = -143.1620228648
γ = 61.960325442</pre>
 
April 2016: This computation took over an hour with MoarVM, running in a VirtualBox linux system guest hosted by a windows laptop with a i7 intel processor.<br>
March 2020: With various improvements to the language, 6.d releases of Raku now run the code 2x to 3x as fast, depending on the hardware used, but even so it is generous to describe it as 'slow'.
 
=={{header|Ruby}}==
Line 1,275 ⟶ 2,598:
Using the standard library Matrix class:
 
<langsyntaxhighlight lang="ruby">require 'matrix'
 
def regression_coefficients y, x
Line 1,282 ⟶ 2,605:
 
(x.t * x).inverse * x.t * y
end</langsyntaxhighlight>
 
Testing 2-dimension:
<langsyntaxhighlight lang="ruby">puts regression_coefficients([1, 2, 3, 4, 5], [ [2, 1, 3, 4, 5] ])</langsyntaxhighlight>
{{out}}
<pre>Matrix[[0.981818181818182]]</pre>
Line 1,291 ⟶ 2,614:
Testing 3-dimension:
Points(x,y,z): [1,1,3], [2,1,4] and [1,2,5]
<langsyntaxhighlight lang="ruby">puts regression_coefficients([3,4,5], [ [1,2,1], [1,1,2] ])</langsyntaxhighlight>
{{out}}
<pre>Matrix[[0.9999999999999996], [2.0]]</pre>
 
=={{header|SPSS}}==
 
First, build a random dataset:
 
<syntaxhighlight lang="text">set rng=mc seed=17760704.
new file.
input program.
vector x(4).
loop #i=1 to 200.
loop #j=1 to 4.
compute x(#j)=rv.normal(0,1).
end loop.
end case.
end loop.
end file.
end input program.
compute y=1.5+0.8*x1-0.7*x2+1.1*x3-1.7*x4+rv.normal(0,1).
execute.</syntaxhighlight>
 
Now use the regression command:
 
<syntaxhighlight lang="text">regression /dependent=y
/method=enter x1 x2 x3 x4.</syntaxhighlight>
 
{{out}}
 
<syntaxhighlight lang="text">Regression
Notes
|--------------------------------------------------------------------|---------------------------------------------------------------------------|
|Output Created |21-MAR-2020 23:17:33 |
|--------------------------------------------------------------------|---------------------------------------------------------------------------|
|Comments | |
|----------------------|---------------------------------------------|---------------------------------------------------------------------------|
|Input |Filter |<none> |
| |---------------------------------------------|---------------------------------------------------------------------------|
| |Weight |<none> |
| |---------------------------------------------|---------------------------------------------------------------------------|
| |Split File |<none> |
| |---------------------------------------------|---------------------------------------------------------------------------|
| |N of Rows in Working Data File |200 |
|----------------------|---------------------------------------------|---------------------------------------------------------------------------|
|Missing Value Handling|Definition of Missing |User-defined missing values are treated as missing. |
| |---------------------------------------------|---------------------------------------------------------------------------|
| |Cases Used |Statistics are based on cases with no missing values for any variable used.|
|--------------------------------------------------------------------|---------------------------------------------------------------------------|
|Syntax |regression /dependent=y /method=enter x1 x2 x3 x4. |
|----------------------|---------------------------------------------|---------------------------------------------------------------------------|
|Resources |Processor Time |00:00:00,00 |
| |---------------------------------------------|---------------------------------------------------------------------------|
| |Elapsed Time |00:00:00,00 |
| |---------------------------------------------|---------------------------------------------------------------------------|
| |Memory Required |4080 bytes |
| |---------------------------------------------|---------------------------------------------------------------------------|
| |Additional Memory Required for Residual Plots|0 bytes |
|------------------------------------------------------------------------------------------------------------------------------------------------|
 
Variables Entered/Removeda
|-----|-----------------|-----------------|------|
|Model|Variables Entered|Variables Removed|Method|
|-----|-----------------|-----------------|------|
|1 |x4, x3, x2, x1b |. |Enter |
|------------------------------------------------|
a Dependent Variable: y
b All requested variables entered.
 
Model Summary
|-----|-----|--------|-----------------|--------------------------|
|Model|R |R Square|Adjusted R Square|Std. Error of the Estimate|
|-----|-----|--------|-----------------|--------------------------|
|1 |,929a|,863 |,860 |,94928 |
|-----------------------------------------------------------------|
a Predictors: (Constant), x4, x3, x2, x1
 
ANOVAa
|----------------|--------------|---|-----------|-------|-----|
|Model |Sum of Squares|df |Mean Square|F |Sig. |
|-----|----------|--------------|---|-----------|-------|-----|
|1 |Regression|1106,659 |4 |276,665 |307,021|,000b|
| |----------|--------------|---|-----------|-------|-----|
| |Residual |175,720 |195|,901 | | |
| |----------|--------------|---|-----------|-------|-----|
| |Total |1282,379 |199| | | |
|-------------------------------------------------------------|
a Dependent Variable: y
b Predictors: (Constant), x4, x3, x2, x1
 
Coefficientsa
|----------------|--------------------------------------|-------------------------|-------|----|
|Model |Unstandardized Coefficients |Standardized Coefficients|t |Sig.|
| |---------------------------|----------|-------------------------| | |
| |B |Std. Error|Beta | | |
|-----|----------|---------------------------|----------|-------------------------|-------|----|
|1 |(Constant)|1,550 |,067 | |23,003 |,000|
| |----------|---------------------------|----------|-------------------------|-------|----|
| |x1 |,831 |,062 |,360 |13,457 |,000|
| |----------|---------------------------|----------|-------------------------|-------|----|
| |x2 |-,604 |,075 |-,215 |-8,051 |,000|
| |----------|---------------------------|----------|-------------------------|-------|----|
| |x3 |1,098 |,065 |,451 |16,989 |,000|
| |----------|---------------------------|----------|-------------------------|-------|----|
| |x4 |-1,770 |,073 |-,656 |-24,306|,000|
|----------------------------------------------------------------------------------------------|
a Dependent Variable: y</syntaxhighlight>
 
=={{header|Stata}}==
Line 1,299 ⟶ 2,726:
First, build a random dataset:
 
<langsyntaxhighlight lang="stata">clear
set seed 17760704
set obs 200
Line 1,305 ⟶ 2,732:
gen x`i'=rnormal()
}
gen y=1.5+0.8*x1-0.7*x2+1.1*x3-1.7*x4+rnormal()</langsyntaxhighlight>
 
Now, use the '''[https://www.stata.com/help.cgi?regress regress]''' command:
 
<syntaxhighlight lang ="stata">reg y x*</langsyntaxhighlight>
 
'''Output'''
Line 1,334 ⟶ 2,761:
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]:
 
<langsyntaxhighlight lang="stata">. di _b[x1]
.75252466
 
Line 1,344 ⟶ 2,771:
 
. di _se[_cons]
.06978623</langsyntaxhighlight>
 
See '''[https://www.stata.com/help.cgi?estatregress_postestimation estat]''',regress '''[https://www.stata.com/help.cgi?predict predictpostestimation]''', '''[https://www.stata.com/help.cgi?estimatesfor estimates]''', '''[https://www.stata.com/help.cgi?margins margins]''' fora exampleslist of commands that can be used after a regression.
 
Here we compute [[wp:Akaike information criterion|Akaike's AIC]], the covariance matrix of the estimates, the predicted values and residuals:
 
<langsyntaxhighlight lang="stata">. estat ic
 
Akaike's information criterion and Bayesian information criterion
Line 1,374 ⟶ 2,801:
 
. predict yhat, xb
. predict r, r</langsyntaxhighlight>
 
=={{header|Tcl}}==
{{tcllib|math::linearalgebra}}
<langsyntaxhighlight lang="tcl">package require math::linearalgebra
namespace eval multipleRegression {
namespace export regressionCoefficients
Line 1,393 ⟶ 2,820:
}
}
namespace import multipleRegression::regressionCoefficients</langsyntaxhighlight>
Using an example from the Wikipedia page on the correlation of height and weight:
<langsyntaxhighlight lang="tcl"># Simple helper just for this example
proc map {n exp list} {
upvar 1 $n v
Line 1,410 ⟶ 2,837:
}
# 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}]]</langsyntaxhighlight>
{{out}} (a 3-vector of coefficients):
<pre>128.81280358170625 -143.16202286630732 61.96032544293041</pre>
 
=={{header|TI-83 BASIC}}==
{{works with|TI-83 BASIC|TI-84Plus 2.55MP}}
<syntaxhighlight lang="ti83b">{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}→L₁
{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₂ </syntaxhighlight>
{{out}}
<pre>
y=ax²+bx+c
a=61.96032544
b=–143.1620229
c=128.8128036</pre>
 
 
=={{header|Ursala}}==
Line 1,418 ⟶ 2,858:
the Lapack library [http://www.netlib.org/lapack/lug/node27.html],
which is callable in Ursala like this:
<langsyntaxhighlight Ursalalang="ursala">regression_coefficients = lapack..dgelsd</langsyntaxhighlight>
test program:
<syntaxhighlight lang="ursala">x =
<lang Ursala>x =
 
<
Line 1,431 ⟶ 2,871:
#cast %eL
 
example = regression_coefficients(x,y)</langsyntaxhighlight>
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 1,445 ⟶ 2,885:
A similar method can be used for regression with complex numbers by substituting
zgelsd for dgelsd, above.
 
=={{header|Visual Basic .NET}}==
{{trans|Java}}
<syntaxhighlight lang="vbnet">Module Module1
 
Sub Swap(Of T)(ByRef x As T, ByRef y As T)
Dim temp = x
x = y
y = temp
End Sub
 
Sub Require(condition As Boolean, message As String)
If condition Then
Return
End If
Throw New ArgumentException(message)
End Sub
 
Class Matrix
Private data As Double(,)
Private rowCount As Integer
Private colCount As Integer
 
Public Sub New(rows As Integer, cols As Integer)
Require(rows > 0, "Need at least one row")
rowCount = rows
 
Require(cols > 0, "Need at least one column")
colCount = cols
 
data = New Double(rows - 1, cols - 1) {}
End Sub
 
Public Sub New(source As Double(,))
Dim rows = source.GetLength(0)
Require(rows > 0, "Need at least one row")
rowCount = rows
 
Dim cols = source.GetLength(1)
Require(cols > 0, "Need at least one column")
colCount = cols
 
data = New Double(rows - 1, cols - 1) {}
For i = 1 To rows
For j = 1 To cols
data(i - 1, j - 1) = source(i - 1, j - 1)
Next
Next
End Sub
 
Default Public Property Index(i As Integer, j As Integer) As Double
Get
Return data(i, j)
End Get
Set(value As Double)
data(i, j) = value
End Set
End Property
 
Public Property Slice(i As Integer) As Double()
Get
Dim m(colCount - 1) As Double
For j = 1 To colCount
m(j - 1) = Index(i, j - 1)
Next
Return m
End Get
Set(value As Double())
Require(colCount = value.Length, "Slice must match the number of columns")
For j = 1 To colCount
Index(i, j - 1) = value(j - 1)
Next
End Set
End Property
 
Public Shared Operator *(m1 As Matrix, m2 As Matrix) As Matrix
Dim rc1 = m1.rowCount
Dim cc1 = m1.colCount
Dim rc2 = m2.rowCount
Dim cc2 = m2.colCount
Require(cc1 = rc2, "Cannot multiply if the first columns does not equal the second rows")
Dim result As New Matrix(rc1, cc2)
For i = 1 To rc1
For j = 1 To cc2
For k = 1 To rc2
result(i - 1, j - 1) += m1(i - 1, k - 1) * m2(k - 1, j - 1)
Next
Next
Next
Return result
End Operator
 
Public Function Transpose() As Matrix
Dim rc = rowCount
Dim cc = colCount
 
Dim trans As New Matrix(cc, rc)
For i = 1 To cc
For j = 1 To rc
trans(i - 1, j - 1) = Index(j - 1, i - 1)
Next
Next
Return trans
End Function
 
Public Sub ToReducedRowEchelonForm()
Dim lead = 0
Dim rc = rowCount
Dim cc = colCount
For r = 1 To rc
If cc <= lead Then
Return
End If
Dim i = r
 
While Index(i - 1, lead) = 0.0
i += 1
If rc = i Then
i = r
lead += 1
If cc = lead Then
Return
End If
End If
End While
 
Dim temp = Slice(i - 1)
Slice(i - 1) = Slice(r - 1)
Slice(r - 1) = temp
 
If Index(r - 1, lead) <> 0.0 Then
Dim div = Index(r - 1, lead)
For j = 1 To cc
Index(r - 1, j - 1) /= div
Next
End If
 
For k = 1 To rc
If k <> r Then
Dim mult = Index(k - 1, lead)
For j = 1 To cc
Index(k - 1, j - 1) -= Index(r - 1, j - 1) * mult
Next
End If
Next
 
lead += 1
Next
End Sub
 
Public Function Inverse() As Matrix
Require(rowCount = colCount, "Not a square matrix")
Dim len = rowCount
Dim aug As New Matrix(len, 2 * len)
For i = 1 To len
For j = 1 To len
aug(i - 1, j - 1) = Index(i - 1, j - 1)
Next
REM augment identity matrix to right
aug(i - 1, i + len - 1) = 1.0
Next
aug.ToReducedRowEchelonForm()
Dim inv As New Matrix(len, len)
For i = 1 To len
For j = len + 1 To 2 * len
inv(i - 1, j - len - 1) = aug(i - 1, j - 1)
Next
Next
Return inv
End Function
End Class
 
Function ConvertArray(source As Double()) As Double(,)
Dim dest(0, source.Length - 1) As Double
For i = 1 To source.Length
dest(0, i - 1) = source(i - 1)
Next
Return dest
End Function
 
Function MultipleRegression(y As Double(), x As Matrix) As Double()
Dim tm As New Matrix(ConvertArray(y))
Dim cy = tm.Transpose
Dim cx = x.Transpose
Return ((x * cx).Inverse * x * cy).Transpose.Slice(0)
End Function
 
Sub Print(v As Double())
Dim it = v.GetEnumerator()
 
Console.Write("[")
If it.MoveNext() Then
Console.Write(it.Current)
End If
While it.MoveNext
Console.Write(", ")
Console.Write(it.Current)
End While
Console.Write("]")
End Sub
 
Sub Main()
Dim y() = {1.0, 2.0, 3.0, 4.0, 5.0}
Dim x As New Matrix({{2.0, 1.0, 3.0, 4.0, 5.0}})
Dim v = MultipleRegression(y, x)
Print(v)
Console.WriteLine()
 
y = {3.0, 4.0, 5.0}
x = New Matrix({
{1.0, 2.0, 1.0},
{1.0, 1.0, 2.0}
})
v = MultipleRegression(y, x)
Print(v)
Console.WriteLine()
 
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}
Dim 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}
 
Dim xs(2, a.Length - 1) As Double
For i = 1 To a.Length
xs(0, i - 1) = 1.0
xs(1, i - 1) = a(i - 1)
xs(2, i - 1) = a(i - 1) * a(i - 1)
Next
x = New Matrix(xs)
v = MultipleRegression(y, x)
Print(v)
Console.WriteLine()
End Sub
 
End Module</syntaxhighlight>
{{out}}
<pre>[0.981818181818182]
[1, 2]
[128.812803579828, -143.162022865304, 61.9603254429854]</pre>
 
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-matrix}}
<syntaxhighlight lang="wren">import "./matrix" for Matrix
 
var multipleRegression = Fn.new { |y, x|
var cy = y.transpose
var cx = x.transpose
return ((x * cx).inverse * x * cy).transpose[0]
}
 
var ys = [
Matrix.new([ [1, 2, 3, 4, 5] ]),
Matrix.new([ [3, 4, 5] ]),
Matrix.new([ [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] ])
]
 
var a = [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]
 
var xs = [
Matrix.new([ [2, 1, 3, 4, 5] ]),
Matrix.new([ [1, 2, 1], [1, 1, 2] ]),
Matrix.new([ List.filled(a.count, 1), a, a.map { |e| e * e }.toList ])
]
 
for (i in 0...ys.count) {
var v = multipleRegression.call(ys[i], xs[i])
System.print(v)
System.print()
}</syntaxhighlight>
 
{{out}}
<pre>
[0.98181818181818]
 
[1, 2]
 
[128.81280357983, -143.1620228653, 61.960325442985]
</pre>
 
=={{header|zkl}}==
Using the GNU Scientific Library:
<langsyntaxhighlight lang="zkl">var [const] GSL=Import("zklGSL"); // libGSL (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 1,456 ⟶ 3,174:
v.format().println();
GSL.Helpers.polyString(v).println();
GSL.Helpers.polyEval(v,height).format().println();</langsyntaxhighlight>
{{out}}
<pre>
Line 1,466 ⟶ 3,184:
Or, using Lists:
{{trans|Common Lisp}}
<langsyntaxhighlight lang="zkl">// Solve a linear system AX=B where A is symmetric and positive definite, so it can be Cholesky decomposed.
fcn linsys(A,B){
n,m:=A.len(),B[1].len(); // A.rows,B.cols
Line 1,521 ⟶ 3,239:
if(M.len()==1) M[0].pump(List,List.create); // 1 row --> n columns
else M[0].zip(M.xplode(1));
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">height:=T(T(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));
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();</langsyntaxhighlight>
{{out}}
<pre>
9,476

edits