Polynomial regression: Difference between revisions

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

 

 

=={{header|Ada}}==

 

function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is

end loop;

return Solve (A * Transpose (A), A * Y);

The function Fit implements least squares approximation of a function defined in the points as specified by the arrays ''x''_''i'' and ''y''_''i''. The basis φ_''j'' is ''x''^''j'', ''j''=0,1,..,''N''. The implementation is straightforward. First the plane matrix A is created. A_ji=φ_''j''(''x''_''i''). Then the linear problem AA^''T''''c''=A''y'' is solved. The result ''c''_''j'' are the coefficients. Constraint_Error is propagated when dimensions of X and Y differ or else when the problem is ill-defined.

===Example===

with Ada.Float_Text_IO; use Ada.Float_Text_IO;

 

Put (C (1), Aft => 3, Exp => 0);

Put (C (2), Aft => 3, Exp => 0);

{{out}}

<pre>

 

<!-- {{does not work with|ELLA ALGOL 68|Any (with appropriate job cards AND formatted transput statements removed) - tested with release 1.8.8d.fc9.i386 - ELLA has no FORMATted transput}} -->

 

MODE

);

print polynomial(d)

{{out}}

<pre>

3x**2+2x+1

1.0848x**2+10.3552x-0.6164

</pre>

 

and fits an order-5 polynomial, so the test data for this task

is hardly challenging!

Max% = 10000

FOR term% = 5 TO 0 STEP -1

PRINT ;vector(term%) " * x^" STR$(term%)

{{out}}

<pre>

 

'''Include''' file (to make the code reusable easily) named <tt>polifitgsl.h</tt>

#define _POLIFITGSL_H

#include <gsl/gsl_multifit.h>

bool polynomialfit(int obs, int degree,

double *dx, double *dy, double *store); /* n, p */

'''Implementation''' (the examples [http://www.gnu.org/software/gsl/manual/html_node/Fitting-Examples.html here] helped alot to code this quickly):

 

bool polynomialfit(int obs, int degree,

return true; /* we do not "analyse" the result (cov matrix mainly)

to know if the fit is "good" */

'''Testing''':

 

#include "polifitgsl.h"

}

return 0;

{{out}}

<pre>1.000000

3.000000</pre>

 

{{libheader|Math.Net}}

{

// Vandermonde matrix

for (int j = 0; j <= degree; j++) v[i, j] = Math.Pow(x[i], j);

var yv = new DenseVector(y).ToColumnMatrix();

// Math.Net doesn't have an "economy" QR, so:

// cut R short to square upper triangle, then recompute Q

var p = r.Inverse().Multiply(q.TransposeThisAndMultiply(yv));

return p.Column(0).ToArray();

Example:

{

const int degree = 2;

Console.WriteLine();

for (int i = 0; i < x.Length; i++ )

Console.ReadKey(true);

 

 

 

=={{header|D}}==

{{trans|Kotlin}}

import std.range;

import std.stdio;

auto y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];

polyRegression(x, y);

{{out}}

<pre>y = 1 + 2x + 3x^2

9 262 262.0

10 321 321.0</pre>

 

=={{header|Emacs Lisp}}==

 

 

 

 

=={{header|Fortran}}==

{{libheader|LAPACK}}

contains

 

end function

 

===Example===

use fitting

implicit none

write (*, '(F9.4)') a

 

 

{{out}} (lower powers first, so this seems the opposite of the Python output):

 

=={{header|FreeBASIC}}==

#macro pivot(num)

#endmacro

For k As Integer=1 To n-1

For row As Integer =k To n-1

For g As Integer=1 To n

Next g

Next row

Next k

For z As Integer=n To 1 Step -1

For j As Integer = n To z+1 Step -1

Next j

 

 

=={{header|GAP}}==

local a;

a := List([0 .. n], i -> List(x, s -> s^i));

# Return coefficients in ascending degree order

PolynomialRegression(x, y, 2);

 

=={{header|gnuplot}}==

 

f(x) = a*x**2 + b*x + c

 

e

 

 

=={{header|Go}}==

===Library gonum/matrix===

 

import (

 

)

 

 

 

 

 

}

 

{{out}}

<pre>

===Library go.matrix===

Least squares solution using QR decomposition and package [http://github.com/skelterjohn/go.matrix go.matrix].

 

import (

}

fmt.Println(c)

{{out}} (lowest order coefficient first)

<pre>

=={{header|Haskell}}==

Uses module Matrix.LU from [http://hackage.haskell.org/package/dsp hackageDB DSP]

import Data.Array

import Control.Monad

polyfit d ry = elems $ solve mat vec where

mat = listArray ((1,1), (d,d)) $ liftM2 concatMap ppoly id [0..fromIntegral $ pred d]

{{out}} in GHCi:

 

=={{header|HicEst}}==

 

x = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)

! called by the solver of the SOLVE function. All variables are global

Theory = p(1)*x(nr)^2 + p(2)*x(nr) + p(3)

{{out}}

<pre>SOLVE performs a (nonlinear) least-square fit (Levenberg-Marquardt):

 

=={{header|Hy}}==

 

(setv x (range 11))

(setv y [1 6 17 34 57 86 121 162 209 262 321])

 

 

=={{header|J}}==

 

(%. ^/~@x:@i.@#) Y

 

Note that this implementation does not use floating point numbers,

but for small problems like this it is inconsequential.

 

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:

(note that this time we used floating point numbers, so that result is approximate rather than exact - it only looks exact because of how J displays floating point numbers (by default, J assumes six digits of accuracy) - changing (i.3) to (x:i.3) would give us an exact result, if that mattered.)

 

=={{header|Julia}}==

The least-squares fit problem for a degree <i>n</i>

can be solved with the built-in backslash operator (coefficients in increasing order of degree):

 

x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]

 

{{out}}

=={{header|Kotlin}}==

{{trans|REXX}}

 

fun polyRegression(x: IntArray, y: IntArray) {

val xm = x.average()

val ym = y.average()

val xym = x.zip(y).map { it.first * it.second }.average()

val x2ym = x.zip(y).map { it.first * it.first * it.second }.average()

println(" Input Approximation")

println(" x y y1")

}

}

 

val x = IntArray(11) { it }

val y = intArrayOf(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)

polyRegression(x, y)

 

{{out}}

</pre>

 

 

Result:

<pre>1 + 2x + 3x^2</pre>

The output of this function is the coefficients of the polynomial which best fit these x,y value pairs.

 

>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];

>> polyfit(x,y,2)

ans =

 

 

=={{header|МК-61/52}}==

Part 1:

ИПC x^2 П2 ИП6 + П6 ИП2 ИПC * ИП7

+ П7 ИП2 x^2 ИП8 + П8 ИПC ИПD *

ИПA + ПA ИП2 ИПD * ИПB + ПB ИПD

 

''Input'': В/О x₁ С/П y₁ С/П x₂ С/П y₂ С/П ...

 

Part 2:

ИПC ИП5 * - ПD ИП4 ИП7 * ИПC ИП6

* - П7 ИП4 ИПA * ИПC ИП9 * -

ИПD ИП8 * ИП7 ИП3 * - / ПB ИПA

ИПB ИП7 * - ИПD / ПA ИП9 ИПB ИП6

 

''Result'': Р9 = a₀, РA = a₁, РB = a₂.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

{{out}}

</pre>

 

 

 

 

 

{{out}}

<pre>

y=3x^2+2x+1

10 321 321

</pre>

 

 

 

 

 

 

=={{header|Python}}==

 

>>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]

>>> coeffs = numpy.polyfit(x,y,deg=2)

>>> coeffs

Substitute back received coefficients.

>>> yf

Find max absolute error:

 

===Example===

For input arrays `x' and `y':

 

>>> print p

2

Thus we confirm once more that for already sorted sequences

the considered quick sort implementation has

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)

 

{{out}}

Alternately, use poly:

 

<pre> (Intercept) poly(x, 2, raw = T)1 poly(x, 2, raw = T)2

1 2 3</pre>

 

=={{header|Racket}}==

#lang racket

(require math plot)

(plot (list (points (map vector xs ys))

(function (poly (fit xs ys 2)))))

{{out}}

[[File:polyreg-racket.png]]

 

=={{header|REXX}}==

* Implementation of http://keisan.casio.com/exec/system/14059932254941

*--------------------------------------------------------------------*/

fun:

Parse Arg x

{{out}}

<pre>y=1+2*x+3*x**2

9 262 262.000

10 321 321.000</pre>

 

=={{header|Ruby}}==

 

def regress x, y, degree

 

mx = Matrix[*x_data]

my = Matrix.column_vector(y)

 

'''Testing:'''

{{out}}

 

=={{header|Sidef}}==

{{trans|Ruby}}

 

 

}

 

2

)

 

{{out}}

 

=={{header|Stata}}==

See '''[http://www.stata.com/help.cgi?fvvarlist Factor variables]''' in Stata help for explanations on the ''c.x##c.x'' syntax.

. input x y

0 1

|

_cons | 1 . . . . .

 

=={{header|Tcl}}==

<!-- This implementation from Emiliano Gavilan;

posted here with his explicit permission -->

 

proc build.matrix {xvec degree} {

set coeffs [math::linearalgebra::solveGauss $A $b]

# show results

This will print:

1.0000000000000207 1.9999999999999958 3.0

which is a close approximation to the correct solution.

 

=={{header|TI-89 BASIC}}==

seq(x,x,0,10) → xs

{1,6,17,34,57,86,121,162,209,262,321} → ys

QuadReg xs,ys

 

<code>seq(''expr'',''var'',''low'',''high'')</code> evaluates ''expr'' with ''var'' bound to integers from ''low'' to ''high'' and returns a list of the results. <code> →</code> is the assignment operator.

whereby the data can be passed as lists rather than arrays,

and all memory management is handled automatically.

#import nat

#import flo

 

test program:

y = <1.,6.,17.,34.,57.,86.,121.,162.,209.,262.,321.>

 

#cast %eL

 

{{out}}

<pre><3.000000e+00,2.000000e+00,1.000000e+00></pre>

 

=={{header|zkl}}==

Using the GNU Scientific Library

xs:=GSL.VectorFromData(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10);

ys:=GSL.VectorFromData(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321);

v.format().println();

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

{{out}}

<pre>

 

Example:

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)

{{out}}<pre>L(1,2,3)</pre>