Polynomial regression

From Rosetta Code
Jump to: navigation, search
Task
Polynomial regression
You are encouraged to solve this task according to the task description, using any language you may know.

Find an approximating polynom of known degree for a given data.

Example: For input data:

x = {0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10};
y = {1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321};

The approximating polynom is:

3 x2 + 2 x + 1

Here, the polynom's coefficients are (3, 2, 1).

This task is intended as a subtask for Measure relative performance of sorting algorithms implementations.


Contents

[edit] Ada

with Ada.Numerics.Real_Arrays;  use Ada.Numerics.Real_Arrays;
 
function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is
A : Real_Matrix (0..N, X'Range); -- The plane
begin
for I in A'Range (2) loop
for J in A'Range (1) loop
A (J, I) := X (I)**J;
end loop;
end loop;
return Solve (A * Transpose (A), A * Y);
end Fit;

The function Fit implements least squares approximation of a function defined in the points as specified by the arrays xi and yi. The basis φj is xj, j=0,1,..,N. The implementation is straightforward. First the plane matrix A is created. Ajij(xi). Then the linear problem AATc=Ay is solved. The result cj are the coefficients. Constraint_Error is propagated when dimensions of X and Y differ or else when the problem is ill-defined.

[edit] Example

with Fit;
with Ada.Float_Text_IO; use Ada.Float_Text_IO;
 
procedure Fitting is
C : constant Real_Vector :=
Fit
( (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0),
(1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0),
2
);
begin
Put (C (0), Aft => 3, Exp => 0);
Put (C (1), Aft => 3, Exp => 0);
Put (C (2), Aft => 3, Exp => 0);
end Fitting;
Output:
 1.000 2.000 3.000

[edit] ALGOL 68

Translation of: Ada
Works with: ALGOL 68 version Standard - lu decomp and lu solve are from the GSL library

Works with: ALGOL 68G version Any - tested with release mk15-0.8b.fc9.i386
MODE FIELD = REAL;
 
MODE
VEC = [0]FIELD,
MAT = [0,0]FIELD;
 
PROC VOID raise index error := VOID: (
print(("stop", new line));
stop
);
 
COMMENT from http://rosettacode.org/wiki/Matrix_Transpose#ALGOL_68 END COMMENT
OP ZIP = ([,]FIELD in)[,]FIELD:(
[2 LWB in:2 UPB in,1 LWB in:1UPB in]FIELD out;
FOR i FROM LWB in TO UPB in DO
out[,i]:=in[i,]
OD;
out
);
 
COMMENT from http://rosettacode.org/wiki/Matrix_multiplication#ALGOL_68 END COMMENT
OP * = (VEC a,b)FIELD: ( # basically the dot product #
FIELD result:=0;
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI;
FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD;
result
);
 
OP * = (VEC a, MAT b)VEC: ( # overload vector times matrix #
[2 LWB b:2 UPB b]FIELD result;
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI;
FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD;
result
);
 
OP * = (MAT a, b)MAT: ( # overload matrix times matrix #
[LWB a:UPB a, 2 LWB b:2 UPB b]FIELD result;
IF 2 LWB a/=LWB b OR 2 UPB a/=UPB b THEN raise index error FI;
FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD;
result
);
 
COMMENT from http://rosettacode.org/wiki/Pyramid_of_numbers#ALGOL_68 END COMMENT
OP / = (VEC a, MAT b)VEC: ( # vector division #
[LWB a:UPB a,1]FIELD transpose a;
transpose a[,1]:=a;
(transpose a/b)[,1]
);
 
OP / = (MAT a, MAT b)MAT:( # matrix division #
[LWB b:UPB b]INT p ;
INT sign;
[,]FIELD lu = lu decomp(b, p, sign);
[LWB a:UPB a, 2 LWB a:2 UPB a]FIELD out;
FOR col FROM 2 LWB a TO 2 UPB a DO
out[,col] := lu solve(b, lu, p, a[,col]) [@LWB out[,col]]
OD;
out
);
 
FORMAT int repr = $g(0)$,
real repr = $g(-7,4)$;
 
PROC fit = (VEC x, y, INT order)VEC:
BEGIN
[0:order, LWB x:UPB x]FIELD a; # the plane #
FOR i FROM 2 LWB a TO 2 UPB a DO
FOR j FROM LWB a TO UPB a DO
a [j, i] := x [i]**j
OD
OD;
( y * ZIP a ) / ( a * ZIP a )
END # fit #;
 
PROC print polynomial = (VEC x)VOID: (
BOOL empty := TRUE;
FOR i FROM UPB x BY -1 TO LWB x DO
IF x[i] NE 0 THEN
IF x[i] > 0 AND NOT empty THEN print ("+") FI;
empty := FALSE;
IF x[i] NE 1 OR i=0 THEN
IF ENTIER x[i] = x[i] THEN
printf((int repr, x[i]))
ELSE
printf((real repr, x[i]))
FI
FI;
CASE i+1 IN
SKIP,print(("x"))
OUT
printf(($"x**"g(0)$,i))
ESAC
FI
OD;
IF empty THEN print("0") FI;
print(new line)
);
 
fitting: BEGIN
VEC c =
fit
( (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0),
(1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0),
2
);
print polynomial(c);
VEC d =
fit
( (0, 1, 2, 3, 4, 5, 6, 7, 8, 9),
(2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0),
2
);
print polynomial(d)
END # fitting #
Output:
3x**2+2x+1
 1.0848x**2+10.3552x-0.6164

[edit] BBC BASIC

The code listed below is good for up to 10000 data points and fits an order-5 polynomial, so the test data for this task is hardly challenging!

      INSTALL @lib$+"ARRAYLIB"
 
Max% = 10000
DIM vector(5), matrix(5,5)
DIM x(Max%), x2(Max%), x3(Max%), x4(Max%), x5(Max%)
DIM x6(Max%), x7(Max%), x8(Max%), x9(Max%), x10(Max%)
DIM y(Max%), xy(Max%), x2y(Max%), x3y(Max%), x4y(Max%), x5y(Max%)
 
npts% = 11
x() = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
y() = 1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321
 
sum_x = SUM(x())
x2() = x() * x()  : sum_x2 = SUM(x2())
x3() = x() * x2()  : sum_x3 = SUM(x3())
x4() = x2() * x2() : sum_x4 = SUM(x4())
x5() = x2() * x3() : sum_x5 = SUM(x5())
x6() = x3() * x3() : sum_x6 = SUM(x6())
x7() = x3() * x4() : sum_x7 = SUM(x7())
x8() = x4() * x4() : sum_x8 = SUM(x8())
x9() = x4() * x5() : sum_x9 = SUM(x9())
x10() = x5() * x5() : sum_x10 = SUM(x10())
 
sum_y = SUM(y())
xy() = x() * y()  : sum_xy = SUM(xy())
x2y() = x2() * y()  : sum_x2y = SUM(x2y())
x3y() = x3() * y()  : sum_x3y = SUM(x3y())
x4y() = x4() * y()  : sum_x4y = SUM(x4y())
x5y() = x5() * y()  : sum_x5y = SUM(x5y())
 
matrix() = \
\ npts%, sum_x, sum_x2, sum_x3, sum_x4, sum_x5, \
\ sum_x, sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, \
\ sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, \
\ sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, \
\ sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, \
\ sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, sum_x10
 
vector() = \
\ sum_y, sum_xy, sum_x2y, sum_x3y, sum_x4y, sum_x5y
 
PROC_invert(matrix())
vector() = matrix().vector()
 
@% = &2040A
PRINT "Polynomial coefficients = "
FOR term% = 5 TO 0 STEP -1
PRINT ;vector(term%) " * x^" STR$(term%)
NEXT
Output:
Polynomial coefficients =
0.0000 * x^5
-0.0000 * x^4
0.0002 * x^3
2.9993 * x^2
2.0012 * x^1
0.9998 * x^0

[edit] C

Include file (to make the code reusable easily) named polifitgsl.h

#ifndef _POLIFITGSL_H
#define _POLIFITGSL_H
#include <gsl/gsl_multifit.h>
#include <stdbool.h>
#include <math.h>
bool polynomialfit(int obs, int degree,
double *dx, double *dy, double *store); /* n, p */
#endif

Implementation (the examples here helped alot to code this quickly):

#include "polifitgsl.h"
 
bool polynomialfit(int obs, int degree,
double *dx, double *dy, double *store) /* n, p */
{
gsl_multifit_linear_workspace *ws;
gsl_matrix *cov, *X;
gsl_vector *y, *c;
double chisq;
 
int i, j;
 
X = gsl_matrix_alloc(obs, degree);
y = gsl_vector_alloc(obs);
c = gsl_vector_alloc(degree);
cov = gsl_matrix_alloc(degree, degree);
 
for(i=0; i < obs; i++) {
gsl_matrix_set(X, i, 0, 1.0);
for(j=0; j < degree; j++) {
gsl_matrix_set(X, i, j, pow(dx[i], j));
}
gsl_vector_set(y, i, dy[i]);
}
 
ws = gsl_multifit_linear_alloc(obs, degree);
gsl_multifit_linear(X, y, c, cov, &chisq, ws);
 
/* store result ... */
for(i=0; i < degree; i++)
{
store[i] = gsl_vector_get(c, i);
}
 
gsl_multifit_linear_free(ws);
gsl_matrix_free(X);
gsl_matrix_free(cov);
gsl_vector_free(y);
gsl_vector_free(c);
return true; /* we do not "analyse" the result (cov matrix mainly)
to know if the fit is "good" */

}

Testing:

#include <stdio.h>
 
#include "polifitgsl.h"
 
#define NP 11
double x[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
double y[] = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321};
 
#define DEGREE 3
double coeff[DEGREE];
 
int main()
{
int i;
 
polynomialfit(NP, DEGREE, x, y, coeff);
for(i=0; i < DEGREE; i++) {
printf("%lf\n", coeff[i]);
}
return 0;
}
Output:
1.000000
2.000000
3.000000

[edit] Common Lisp

Uses the routine (lsqr A b) from Multiple regression and (mtp A) from Matrix transposition.

;; Least square fit of a polynomial of order n the x-y-curve.
(defun polyfit (x y n)
(let* ((m (cadr (array-dimensions x)))
(A (make-array `(,m ,(+ n 1)) :initial-element 0)))
(loop for i from 0 to (- m 1) do
(loop for j from 0 to n do
(setf (aref A i j)
(expt (aref x 0 i) j))))
(lsqr A (mtp y))))

Example:

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

[edit] C#

Library: Math.Net
        public static double[] Polyfit(double[] x, double[] y, int degree)
{
// Vandermonde matrix
var v = new DenseMatrix(x.Length, degree + 1);
for (int i = 0; i < v.RowCount; i++)
for (int j = 0; j <= degree; j++) v[i, j] = Math.Pow(x[i], j);
var yv = new DenseVector(y).ToColumnMatrix();
QR qr = v.QR();
// Math.Net doesn't have an "economy" QR, so:
// cut R short to square upper triangle, then recompute Q
var r = qr.R.SubMatrix(0, degree + 1, 0, degree + 1);
var q = v.Multiply(r.Inverse());
var p = r.Inverse().Multiply(q.TransposeThisAndMultiply(yv));
return p.Column(0).ToArray();
}

Example:

        static void Main(string[] args)
{
const int degree = 2;
var x = new[] {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
var y = new[] {1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0};
var p = Polyfit(x, y, degree);
foreach (var d in p) Console.Write("{0} ",d);
Console.WriteLine();
for (int i = 0; i < x.Length; i++ )
Console.WriteLine("{0} => {1} diff {2}", x[i], Polyval(p,x[i]), y[i] - Polyval(p,x[i]));
Console.ReadKey(true);
}

[edit] Fortran

Library: LAPACK
module fitting
contains
 
function polyfit(vx, vy, d)
implicit none
integer, intent(in) :: d
integer, parameter :: dp = selected_real_kind(15, 307)
real(dp), dimension(d+1) :: polyfit
real(dp), dimension(:), intent(in) :: vx, vy
 
real(dp), dimension(:,:), allocatable :: X
real(dp), dimension(:,:), allocatable :: XT
real(dp), dimension(:,:), allocatable :: XTX
 
integer :: i, j
 
integer :: n, lda, lwork
integer :: info
integer, dimension(:), allocatable :: ipiv
real(dp), dimension(:), allocatable :: work
 
n = d+1
lda = n
lwork = n
 
allocate(ipiv(n))
allocate(work(lwork))
allocate(XT(n, size(vx)))
allocate(X(size(vx), n))
allocate(XTX(n, n))
 
! prepare the matrix
do i = 0, d
do j = 1, size(vx)
X(j, i+1) = vx(j)**i
end do
end do
 
XT = transpose(X)
XTX = matmul(XT, X)
 
! calls to LAPACK subs DGETRF and DGETRI
call DGETRF(n, n, XTX, lda, ipiv, info)
if ( info /= 0 ) then
print *, "problem"
return
end if
call DGETRI(n, XTX, lda, ipiv, work, lwork, info)
if ( info /= 0 ) then
print *, "problem"
return
end if
 
polyfit = matmul( matmul(XTX, XT), vy)
 
deallocate(ipiv)
deallocate(work)
deallocate(X)
deallocate(XT)
deallocate(XTX)
 
end function
 
end module

[edit] Example

program PolynomalFitting
use fitting
implicit none
 
! let us test it
integer, parameter :: degree = 2
integer, parameter :: dp = selected_real_kind(15, 307)
integer :: i
real(dp), dimension(11) :: x = (/ (i,i=0,10) /)
real(dp), dimension(11) :: y = (/ 1, 6, 17, 34, &
57, 86, 121, 162, &
209, 262, 321 /)
real(dp), dimension(degree+1) :: a
 
a = polyfit(x, y, degree)
 
write (*, '(F9.4)') a
 
end program
Output:
(lower powers first, so this seems the opposite of the Python output):
   1.0000
   2.0000
   3.0000

[edit] freebasic

 
Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)
Dim As Integer n=Ubound(matrix,1)
Redim ans(0):Redim ans(1 To n)
Dim As Double b(1 To n,1 To n),r(1 To n)
For c As Integer=1 To n 'take copies
r(c)=rhs(c)
For d As Integer=1 To n
b(c,d)=matrix(c,d)
Next d
Next c
#macro pivot(num)
For p1 As Integer = num To n - 1
For p2 As Integer = p1 + 1 To n
If Abs(b(p1,num))<Abs(b(p2,num)) Then
Swap r(p1),r(p2)
For g As Integer=1 To n
Swap b(p1,g),b(p2,g)
Next g
End If
Next p2
Next p1
#endmacro
For k As Integer=1 To n-1
pivot(k) 'full pivoting
For row As Integer =k To n-1
If b(row+1,k)=0 Then Exit For
Var f=b(k,k)/b(row+1,k)
r(row+1)=r(row+1)*f-r(k)
For g As Integer=1 To n
b((row+1),g)=b((row+1),g)*f-b(k,g)
Next g
Next row
Next k
'back substitute
For z As Integer=n To 1 Step -1
ans(z)=r(z)/b(z,z)
For j As Integer = n To z+1 Step -1
ans(z)=ans(z)-(b(z,j)*ans(j)/b(z,z))
Next j
Next z
End Sub
 
'Interpolate through points.
Sub Interpolate(x_values() As Double,y_values() As Double,p() As Double)
Var n=Ubound(x_values)
Redim p(0):Redim p(1 To n)
Dim As Double matrix(1 To n,1 To n),rhs(1 To n)
For a As Integer=1 To n
rhs(a)=y_values(a)
For b As Integer=1 To n
matrix(a,b)=x_values(a)^(b-1)
Next b
Next a
'Solve the linear equations
GaussJordan(matrix(),rhs(),p())
End Sub
 
'======================== SET UP THE POINTS ===============
 
Dim As Double x(1 To ...)={0,1,2,3,4,5,6,7,8,9,10}
Dim As Double y(1 To ...)={1,6,17,34,57,86,121,162,209,262,321}
 
Redim As Double Poly(0)
'Get the polynomial Poly()
Interpolate(x(),y(),Poly())
 
'print coefficients to console
print "Polynomial Coefficients:"
print
For z As Integer=1 To Ubound(Poly)
If z=1 Then
Print "constant term ";tab(20);Poly(z)
Else
Print tab(8); "x^";z-1;" = ";tab(20);Poly(z)
End If
Next z
 
sleep
 
Output:
Polynomial Coefficients:

constant term       1
       x^ 1 =       2
       x^ 2 =       3
       x^ 3 =       0
       x^ 4 =       0
       x^ 5 =       0
       x^ 6 =       0
       x^ 7 =       0
       x^ 8 =       0
       x^ 9 =       0
       x^ 10 =      0


[edit] GAP

PolynomialRegression := function(x, y, n)
local a;
a := List([0 .. n], i -> List(x, s -> s^i));
return TransposedMat((a * TransposedMat(a))^-1 * a * TransposedMat([y]))[1];
end;
 
x := [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
y := [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
 
# Return coefficients in ascending degree order
PolynomialRegression(x, y, 2);
# [ 1, 2, 3 ]

[edit] gnuplot

# The polynomial approximation
f(x) = a*x**2 + b*x + c
 
# Initial values for parameters
a = 0.1
b = 0.1
c = 0.1
 
# Fit f to the following data by modifying the variables a, b, c
fit f(x) '-' via a, b, c
0 1
1 6
2 17
3 34
4 57
5 86
6 121
7 162
8 209
9 262
10 321
e
 
print sprintf("\n --- \n Polynomial fit: %.4f x^2 + %.4f x + %.4f\n", a, b, c)

[edit] Go

Least squares solution using QR decomposition and package gomatrix.

package main
 
import (
"code.google.com/p/gomatrix/matrix"
"fmt"
)
 
var xGiven = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
var yGiven = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
var degree = 2
 
func main() {
m := len(yGiven)
n := degree + 1
y := matrix.MakeDenseMatrix(yGiven, m, 1)
x := matrix.Zeros(m, n)
for i := 0; i < m; i++ {
ip := float64(1)
for j := 0; j < n; j++ {
x.Set(i, j, ip)
ip *= xGiven[i]
}
}
 
q, r := x.QR()
qty, err := q.Transpose().Times(y)
if err != nil {
fmt.Println(err)
return
}
c := make([]float64, n)
for i := n - 1; i >= 0; i-- {
c[i] = qty.Get(i, 0)
for j := i + 1; j < n; j++ {
c[i] -= c[j] * r.Get(i, j)
}
c[i] /= r.Get(i, i)
}
fmt.Println(c)
}
Output:
(lowest order coefficient first)
[0.9999999999999758 2.000000000000015 2.999999999999999]

[edit] Haskell

Uses module Matrix.LU from hackageDB DSP

import Data.List
import Data.Array
import Control.Monad
import Control.Arrow
import Matrix.LU
 
ppoly p x = map (x**) p
 
polyfit d ry = elems $ solve mat vec where
mat = listArray ((1,1), (d,d)) $ liftM2 concatMap ppoly id [0..fromIntegral $ pred d]
vec = listArray (1,d) $ take d ry
Output:
in GHCi:
*Main> polyfit 3 [1,6,17,34,57,86,121,162,209,262,321]
[1.0,2.0,3.0]

[edit] HicEst

REAL :: n=10, x(n), y(n), m=3, p(m)
 
x = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y = (1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
 
p = 2 ! initial guess for the polynom's coefficients
 
SOLVE(NUL=Theory()-y(nr), Unknown=p, DataIdx=nr, Iters=iterations)
 
WRITE(ClipBoard, Name) p, iterations
 
FUNCTION Theory()
! called by the solver of the SOLVE function. All variables are global
Theory = p(1)*x(nr)^2 + p(2)*x(nr) + p(3)
END
Output:
SOLVE performs a (nonlinear) least-square fit (Levenberg-Marquardt): 
p(1)=2.997135145; p(2)=2.011348347; p(3)=0.9906627242; iterations=19;

[edit] Hy

(import [numpy [polyfit]])
 
(setv x (range 11))
(setv y [1 6 17 34 57 86 121 162 209 262 321])
 
(print (polyfit x y 2))

[edit] J

   X=:i.# Y=:1 6 17 34 57 86 121 162 209 262 321
Y (%. (^/ x:@i.@#)) X
1 2 3 0 0 0 0 0 0 0 0

Note that this implementation does not use floating point numbers, so we do not introduce floating point errors. Using exact arithmetic has a speed penalty, but for small problems like this it is inconsequential.

The above solution fits a polynomial of order 11. If the order of the polynomial is known to be 3 (as is implied in the task description) then the following solution is probably preferable:

   Y (%. (i.3) ^/~ ]) X
1 2 3

[edit] Julia

The least-squares fit problem for a degree n can be solved with the built-in backslash operator:

function polyfit(x, y, n)
A = [ float(x[i])^p for i = 1:length(x), p = 0:n ]
A \ y
end
Output:
julia> x = [0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10]
julia> y = [1,  6,  17, 34, 57, 86, 121, 162, 209, 262, 321]
julia> polyfit(x, y, 2)
3-element Array{Float64,1}:
 1.0
 2.0
 3.0

(giving the coefficients in increasing order of degree).

[edit] Mathematica

Using the built-in "Fit" function.

data = Transpose@{Range[0, 10], {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}};
Fit[data, {1, x, x^2}, x]

Result:

1 + 2x + 3x^2

[edit] MATLAB

Matlab has a built-in function "polyfit(x,y,n)" which performs this task. The arguments x and y are vectors which are parametrized by the index suck that pointi = (xi,yi) and the argument n is the order of the polynomial you want to fit. The output of this function is the coefficients of the polynomial which best fit these x,y value pairs.

>> x = [0,  1,  2,  3,  4,  5,  6,   7,   8,   9,   10];
>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
>> polyfit(x,y,2)
 
ans =
 
2.999999999999998 2.000000000000019 0.999999999999956

[edit] МК-61/52

Part 1:

ПC	С/П	ПD	ИП9	+	П9	ИПC	ИП5	+	П5
ИПC x^2 П2 ИП6 + П6 ИП2 ИПC * ИП7
+ П7 ИП2 x^2 ИП8 + П8 ИПC ИПD *
ИПA + ПA ИП2 ИПD * ИПB + ПB ИПD
КИП4 С/П БП 00

Input: В/О x1 С/П y1 С/П x2 С/П y2 С/П ...

Part 2:

ИП5	ПC	ИП6	ПD	П2	ИП7	П3	ИП4	ИПD	*
ИПC ИП5 * - ПD ИП4 ИП7 * ИПC ИП6
* - П7 ИП4 ИПA * ИПC ИП9 * -
ПA ИП4 ИП3 * ИП2 ИП5 * - П3 ИП4
ИП8 * ИП2 ИП6 * - П8 ИП4 ИПB *
ИП2 ИП9 * - ИПD * ИП3 ИПA * -
ИПD ИП8 * ИП7 ИП3 * - / ПB ИПA
ИПB ИП7 * - ИПD / ПA ИП9 ИПB ИП6
* - ИПA ИП5 * - ИП4 / П9 С/П

Result: Р9 = a0, РA = a1, РB = a2.

[edit] Octave

x = [0:10];
y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
coeffs = polyfit(x, y, 2)

[edit] PARI/GP

Lagrange interpolating polynomial:

polinterpolate([0,1,2,3,4,5,6,7,8,9,10],[1,6,17,34,57,86,121,162,209,262,321])
Output:
3*x^2 + 2*x + 1

Least-squares fit:

V=[1,6,17,34,57,86,121,162,209,262,321]~;
M=matrix(#V,3,i,j,(i-1)^(j-1));Polrev(matsolve(M~*M,M~*V))
Output:
3*x^2 + 2*x + 1

Code thanks to Bill Allombert

[edit] Python

Library: numpy
>>> x = [0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10]
>>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]
>>> coeffs = numpy.polyfit(x,y,deg=2)
>>> coeffs
array([ 3., 2., 1.])

Substitute back received coefficients.

>>> yf = numpy.polyval(numpy.poly1d(coeffs), x)
>>> yf
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.])

Find max absolute error:

>>> '%.1g' % max(y-yf)
'1e-013'

[edit] Example

For input arrays `x' and `y':

>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0]
>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N')
>>> print p
2
1.085 N + 10.36 N - 0.6164

Thus we confirm once more that for already sorted sequences the considered quick sort implementation has quadratic dependence on sequence length (see Example section for Python language on Query Performance page).

[edit] R

The easiest (and most robust) approach to solve this in R is to use the base package's lm function which will find the least squares solution via a QR decomposition:

 
x <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <- c(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
coef(lm(y ~ x + I(x^2)))
Output:
(Intercept)           x      I(x^2) 
          1           2           3 

[edit] Racket

 
#lang racket
(require math plot)
 
(define xs '(0 1 2 3 4 5 6 7 8 9 10))
(define ys '(1 6 17 34 57 86 121 162 209 262 321))
 
(define (fit x y n)
(define Y (->col-matrix y))
(define V (vandermonde-matrix x (+ n 1)))
(define VT (matrix-transpose V))
(matrix->vector (matrix-solve (matrix* VT V) (matrix* VT Y))))
 
(define ((poly v) x)
(for/sum ([c v] [i (in-naturals)])
(* c (expt x i))))
 
(plot (list (points (map vector xs ys))
(function (poly (fit xs ys 2)))))
 
Output:

Polyreg-racket.png

[edit] Ruby

require 'matrix'
 
def regress x, y, degree
x_data = x.map { |xi| (0..degree).map { |pow| (xi**pow).to_f } }
 
mx = Matrix[*x_data]
my = Matrix.column_vector(y)
 
((mx.t * mx).inv * mx.t * my).transpose.to_a[0]
end

Testing:

betas = regress [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321],
2
 
p betas
Output:
[1.00000000000018, 2.00000000000011, 3.00000000000001]

[edit] Tcl

Library: Tcllib (Package: math::linearalgebra)
package require math::linearalgebra
 
proc build.matrix {xvec degree} {
set sums [llength $xvec]
for {set i 1} {$i <= 2*$degree} {incr i} {
set sum 0
foreach x $xvec {
set sum [expr {$sum + pow($x,$i)}]
}
lappend sums $sum
}
 
set order [expr {$degree + 1}]
set A [math::linearalgebra::mkMatrix $order $order 0]
for {set i 0} {$i <= $degree} {incr i} {
set A [math::linearalgebra::setrow A $i [lrange $sums $i $i+$degree]]
}
return $A
}
 
proc build.vector {xvec yvec degree} {
set sums [list]
for {set i 0} {$i <= $degree} {incr i} {
set sum 0
foreach x $xvec y $yvec {
set sum [expr {$sum + $y * pow($x,$i)}]
}
lappend sums $sum
}
 
set x [math::linearalgebra::mkVector [expr {$degree + 1}] 0]
for {set i 0} {$i <= $degree} {incr i} {
set x [math::linearalgebra::setelem x $i [lindex $sums $i]]
}
return $x
}
 
# Now, to solve the example from the top of this page
set x {0 1 2 3 4 5 6 7 8 9 10}
set y {1 6 17 34 57 86 121 162 209 262 321}
 
# build the system A.x=b
set degree 2
set A [build.matrix $x $degree]
set b [build.vector $x $y $degree]
# solve it
set coeffs [math::linearalgebra::solveGauss $A $b]
# show results
puts $coeffs

This will print:

1.0000000000000207 1.9999999999999958 3.0

which is a close approximation to the correct solution.

[edit] TI-89 BASIC

DelVar x
seq(x,x,0,10) → xs
{1,6,17,34,57,86,121,162,209,262,321} → ys
QuadReg xs,ys
Disp regeq(x)

seq(expr,var,low,high) evaluates expr with var bound to integers from low to high and returns a list of the results. is the assignment operator. QuadReg, "quadratic regression", does the fit and stores the details in a number of standard variables, including regeq, which receives the fitted quadratic (polynomial) function itself. We then apply that function to the (undefined as ensured by DelVar) variable x to obtain the expression in terms of x, and display it.

Output:

3.·x2 + 2.·x + 1.

[edit] Ursala

Library: LAPACK

The fit function defined below returns the coefficients of an nth-degree polynomial in order of descending degree fitting the lists of inputs x and outputs y. The real work is done by the dgelsd function from the lapack library. Ursala provides a simplified interface to this library whereby the data can be passed as lists rather than arrays, and all memory management is handled automatically.

#import std
#import nat
#import flo
 
(fit "n") ("x","y") = ..dgelsd\"y" (gang \/*pow float*x iota successor "n")* "x"

test program:

x = <0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.>
y = <1.,6.,17.,34.,57.,86.,121.,162.,209.,262.,321.>
 
#cast %eL
 
example = fit2(x,y)
Output:
<3.000000e+00,2.000000e+00,1.000000e+00>

[edit] zkl

Translation of: Common Lisp

Uses the code from Multiple regression#zkl.

Example:

polyfit(T(T(0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0)), 
T(T(1.0,6.0,17.0,34.0,57.0,86.0,121.0,162.0,209.0,262.0,321.0)), 2)
.flatten().println();
Output:
L(1,2,3)
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox