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.
Ada
<lang 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;</lang> 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. Aji=φj(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.
Example
<lang ada>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;</lang> Sample output:
1.000 2.000 3.000
ALGOL 68
<lang algol68>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 #</lang> Output:
3x**2+2x+1 1.0848x**2+10.3552x-0.6164
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! <lang bbcbasic> 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</lang>
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
C
Include file (to make the code reusable easily) named polifitgsl.h <lang c>#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</lang>
Implementation (the examples here helped alot to code this quickly): <lang c>#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" */ }</lang> Testing: <lang c>#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;
}</lang> Output of the test:
1.000000 2.000000 3.000000
Common Lisp
Uses the routine (lsqr A b) from Multiple regression and (mtp A) from Matrix transposition.
<lang lisp>;; 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))))</lang>
Example:
<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))</lang>
C#
<lang C sharp> 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(); }</lang>
Example: <lang C sharp> 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); }</lang>
Fortran
<lang fortran>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</lang>
Example
<lang fortran>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</lang>
Output (lower powers first, so this seems the opposite of the Python output):
1.0000 2.0000 3.0000
GAP
<lang 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 ]</lang>
gnuplot
<lang 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)</lang>
Go
Least squares solution using QR decomposition and package gomatrix. <lang go>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)
}</lang> Output (lowest order coefficient first)
[0.9999999999999758 2.000000000000015 2.999999999999999]
Haskell
Uses module Matrix.LU from hackageDB DSP <lang haskell>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</lang>
Output in GHCi: <lang haskell>*Main> polyfit 3 [1,6,17,34,57,86,121,162,209,262,321] [1.0,2.0,3.0]</lang>
HicEst
<lang 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</lang>
<lang hicest>SOLVE performs a (nonlinear) least-square fit (Levenberg-Marquardt): p(1)=2.997135145; p(2)=2.011348347; p(3)=0.9906627242; iterations=19;</lang>
J
<lang 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</lang>
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. <lang j> Y (%. (i.3) ^/~ ]) X 1 2 3</lang>
Mathematica
Using the built-in "Fit" function.
<lang Mathematica>data = Transpose@{Range[0, 10], {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}}; Fit[data, {1, x, x^2}, x]</lang> Result:
1 + 2x + 3x^2
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 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.
<lang MATLAB>>> 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</lang>
Octave
<lang octave>x = [0:10]; y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]; coeffs = polyfit(x, y, 2)</lang>
PARI/GP
Lagrange interpolating polynomial: <lang parigp>polinterpolate([0,1,2,3,4,5,6,7,8,9,10],[1,6,17,34,57,86,121,162,209,262,321])</lang> Output:
3*x^2 + 2*x + 1
Least-squares fit: <lang parigp>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))</lang> Output:
3*x^2 + 2*x + 1
Code thanks to Bill Allombert
Python
<lang python>>>> 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.])</lang> Substitute back received coefficients. <lang python>>>> yf = numpy.polyval(numpy.poly1d(coeffs), x) >>> yf array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.])</lang> Find max absolute error. <lang python>>>> '%.1g' % max(y-yf) '1e-013'</lang>
Example
For input arrays `x' and `y': <lang python>>>> 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]</lang>
<lang python>>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N') >>> print p
2
1.085 N + 10.36 N - 0.6164</lang> 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).
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:
<lang R> 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)))</lang>
Output <lang R> (Intercept) x I(x^2)
1 2 3
</lang>
Racket
<lang 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)))))
</lang> Output: <lang racket> </lang>
Ruby
<lang 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</lang> Testing: <lang ruby>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</lang> Output:
[1.00000000000018, 2.00000000000011, 3.00000000000001]
Tcl
<lang tcl>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</lang> This will print:
1.0000000000000207 1.9999999999999958 3.0
which is a close approximation to the correct solution.
TI-89 BASIC
<lang ti89b>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)</lang>
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.
Ursala
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. <lang Ursala>#import std
- import nat
- import flo
(fit "n") ("x","y") = ..dgelsd\"y" (gang \/*pow float*x iota successor "n")* "x"</lang> test program: <lang Ursala>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)</lang> output:
<3.000000e+00,2.000000e+00,1.000000e+00>