Multiple regression: Difference between revisions

 

matrices.ads:

type Element_Type is private;

Zero : Element_Type;

Not_Square_Matrix : exception;

Not_Invertible : exception;

 

matrices.adb:

function "*" (Left, Right : Matrix) return Matrix is

Result : Matrix (Left'Range (1), Right'Range (2)) :=

return Result;

end Transpose;

 

Example multiple_regression.adb:

with Matrices;

procedure Multiple_Regression is

Output_Matrix (Float_Matrices.To_Matrix (Coefficients));

end;

 

{{out}}

=={{header|BBC BASIC}}==

{{works with|BBC BASIC for Windows}}

INSTALL @lib$+"ARRAYLIB"

t() = t().m()

c() = y().t()

{{out}}

<pre>

 

=={{header|C}}==

#include <gsl/gsl_matrix.h>

#include <gsl/gsl_math.h>

gsl_multifit_linear_free(wspc);

 

 

=={{header|C++}}==

{{trans|Java}}

#include <iostream>

 

 

return 0;

{{out}}

<pre>[0.981818]

=={{header|C sharp|C#}}==

{{libheader|Math.Net}}

using MathNet.Numerics.LinearRegression;

using MathNet.Numerics.LinearAlgebra;

Console.WriteLine(β);

}

 

{{out}}

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)

(linsys (mmul (mtp A) A)

(mmul (mtp A) b)))

 

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))

 

=={{header|D}}==

{{trans|Java}}

import std.array;

import std.exception;

v = multipleRegression(y, x);

v.writeln;

{{out}}

<pre>[0.981818]

{{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

 

{{out}}

 

=={{header|ERRE}}==

 

!$DOUBLE

END FOR

 

{{out}}

<pre>LINEAR SYSTEM COEFFICENTS

{{libheader|SLATEC}} [http://netlib.org/slatec/ Available at the Netlib]

 

* MR - multiple regression using the SLATEC library routine DHFTI

*

STOP 'program complete'

END

{{out}}

<pre>

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 (

x, y := givens()

fmt.Printf("%.4f\n", mat64.Formatted(mat64.QR(x).Solve(y)))

{{out}}

<pre>

</pre>

===Library go.matrix===

 

import (

}

fmt.Println(c)

{{out}}

<pre>

=={{header|Haskell}}==

Using package [http://hackage.haskell.org/package/hmatrix hmatrix] from HackageDB

import Numeric.LinearAlgebra.LAPACK

 

v :: Matrix Double

v = (3><1)

Using lapack::dgels

(3><1)

[ 0.9335611922087276

, 1.101323491272865

Or

(3><1)

[ 0.9335611922087278

, 1.101323491272865

 

=={{header|Hy}}==

[numpy [ones column-stack]]

[numpy.random [randn]]

(print (first (lstsq

(column-stack (, (ones n) x1 x2 (* x1 x2)))

 

=={{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

 

Breaking it down:

X=: (x^0) ,. (x^1) ,. (x^2) NB. equivalent of previous line

4{.X NB. show first 4 rows of X

NB. y %. X does matrix division and gives the regression coefficients

y %. X

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

X (%.@:xpy@[ mp xpy) y

128.814 _143.163 61.9606

 

LAPACK routines are also available via the Addon <tt>math/lapack</tt>.

load 'math/lapack/gels'

gels_jlapack_ X;y

 

=={{header|Java}}==

{{trans|Kotlin}}

import java.util.Objects;

 

printVector(v);

}

{{out}}

<pre>[0.9818181818181818]

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) {

)

);

{{out}}

<pre>0.9818181818181818

 

'''Preliminaries'''

reduce range(0;a|length) as $i (0; . + (a[$i] * b[$i]) );

 

reduce range(0; $p) as $j

(.;

 

'''Multiple Regression'''

(y|transpose) as $cy

| (x|transpose) as $cx

 

range(0; ys|length) as $i

{{out}}

<pre>

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];

{{out}}

<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.

 

typealias Vector = DoubleArray

v = multipleRegression(y, x)

printVector(v)

 

{{out}}

First build a random dataset.

 

X:=<ArrayTools[RandomArray](n,4,distribution=normal)|Vector(n,1,datatype=float[8])>:

 

Now the linear regression, with either the LinearAlgebra package, or the Statistics package.

 

# [4.33701132468683, -2.78654498997457, -1.41840666085642, 2.92065133466547, 1.76076442997642]

 

# R-squared: 0.9767, Adjusted R-squared: 0.9761

# 4.33701132468683 x1 - 2.78654498997457 x2 - 1.41840666085642 x3

 

=={{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};

{{out}}

<pre>{128.813, -143.162, 61.9603}</pre>

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

 

In its transposed form yt = Xt * bt, the backslash operator can be used.

 

bt = Xt \ yt

bt =

-0.0561589

-0.1752146

 

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

 

 

Instead of "/", the slash operator, one can also write :

or

 

=={{header|Nim}}==

{{libheader|arraymancer}}

 

import math

var a = stack(height.ones_like, height, height *. height, axis = 1)

 

 

{{out}}

 

=={{header|PARI/GP}}==

addhelp(pseudoinv, "pseudoinv(M): Moore pseudoinverse of the matrix M.");

 

 

=={{header|Perl}}==

use warnings;

use Statistics::Regression;

my @coeff = $reg->theta();

 

{{out}}

<pre>const 128.813

=={{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>

<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>

 

=={{header|PicoLisp}}==

 

# Matrix transposition

(car X) ) ) ) )

(T (> (inc 'Lead) Cols)) ) )

{{trans|JavaScript}}

(let N (length Mat)

(unless (= N (length (car Mat)))

X (columnVector (2.0 1.0 3.0 4.0 5.0)) )

 

{{out}}

<pre>-> "0.98181818181818182"</pre>

{{libheader|NumPy}}

'''Method with matrix operations'''

 

height = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,

y = np.mat(weight)

 

{{out}}

<pre>

</pre>

'''Using numpy lstsq function'''

 

height = [1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63,

y = weight

 

{{out}}

<pre>

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)

 

 

{{out}}

is useful to illustrate R's model description and linear algebra capabilities.

 

 

## parse and evaluate the model formula

}

 

 

This produces the same coefficients as lm()

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)

 

# (Intercept) x I(x^2)

 

=={{header|Racket}}==

#lang racket

(require math)

(define (fit X y)

(matrix-solve (matrix* (T X) X) (matrix* (T X) y)))

Test:

(fit (matrix [[1 2]

[2 5]

{{out}}

(array #[#[9 1/3] #[-3 1/3]])

 

=={{header|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>;

say "α = ", ($w∧$h1∧$h2)·$I.reversion/$I2;

say "β = ", ($w∧$h2∧$h0)·$I.reversion/$I2;

{{out}}

<pre>α = 128.81280357844

Using the standard library Matrix class:

 

 

def regression_coefficients y, x

 

(x.t * x).inverse * x.t * y

 

Testing 2-dimension:

{{out}}

<pre>Matrix[[0.981818181818182]]</pre>

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>

First, build a random dataset:

 

new file.

input program.

end input program.

compute y=1.5+0.8*x1-0.7*x2+1.1*x3-1.7*x4+rv.normal(0,1).

 

Now use the regression command:

 

 

{{out}}

 

Notes

|--------------------------------------------------------------------|---------------------------------------------------------------------------|

| |x4 |-1,770 |,073 |-,656 |-24,306|,000|

|----------------------------------------------------------------------------------------------|

 

=={{header|Stata}}==

First, build a random dataset:

 

set seed 17760704

set obs 200

gen x`i'=rnormal()

}

 

Now, use the '''[https://www.stata.com/help.cgi?regress regress]''' command:

 

 

'''Output'''

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

 

 

. di _se[_cons]

 

See '''[https://www.stata.com/help.cgi?regress_postestimation regress postestimation]''' for a list 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:

 

 

Akaike's information criterion and Bayesian information criterion

 

. predict yhat, xb

 

=={{header|Tcl}}==

{{tcllib|math::linearalgebra}}

namespace eval multipleRegression {

namespace export regressionCoefficients

}

}

Using an example from the Wikipedia page on the correlation of height and weight:

proc map {n exp list} {

upvar 1 $n v

}

# Wikipedia states that fitting up to the square of x[i] is worth it

{{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}}

{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₂

{{out}}

<pre>

the Lapack library [http://www.netlib.org/lapack/lug/node27.html],

which is callable in Ursala like this:

test program:

 

<

#cast %eL

 

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,

=={{header|Visual Basic .NET}}==

{{trans|Java}}

 

Sub Swap(Of T)(ByRef x As T, ByRef y As T)

End Sub

 

{{out}}

<pre>[0.981818181818182]

{{trans|Kotlin}}

{{libheader|Wren-matrix}}

 

var multipleRegression = Fn.new { |y, x|

System.print(v)

System.print()

 

{{out}}

=={{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);

v.format().println();

GSL.Helpers.polyString(v).println();

{{out}}

<pre>

Or, using Lists:

{{trans|Common Lisp}}

fcn linsys(A,B){

n,m:=A.len(),B[1].len(); // A.rows,B.cols

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));

{{out}}

<pre>