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

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

{{libheader|Math.Net}}

{

// Vandermonde matrix

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

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

Example:

{

const int degree = 2;

Console.WriteLine("{0} => {1} diff {2}", x[i], Polynomial.Evaluate(x[i], p), y[i] - Polynomial.Evaluate(x[i], p));

Console.ReadKey(true);

 

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

{{trans|Java}}

#include <iostream>

#include <numeric>

 

return 0;

{{out}}

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

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

 

(defun polyfit (x y n)

(let* ((m (cadr (array-dimensions x)))

(setf (aref A i j)

(expt (aref x 0 i) j))))

 

Example:

 

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

 

{{trans|D}}

{{works with|Java|8}}

import java.util.function.IntToDoubleFunction;

import java.util.stream.IntStream;

private static void polyRegression(int[] x, int[] y) {

int n = x.length;

double xm = Arrays.stream(x).average().orElse(Double.NaN);

double ym = Arrays.stream(y).average().orElse(Double.NaN);

double xym = 0.0;

for (int i = 0; i < x.length && i < y.length; ++i) {

polyRegression(x, y);

}

{{out}}

<pre>y = 1.0 + 2.0x + 3.0x^2

9 262 262.0

10 321 321.0</pre>

 

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

=={{header|Lua}}==

{{trans|Modula-2}}

return a + (b + c * x) * x

end

local xa = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}

local ya = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}

{{out}}

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

 

=={{header|Maple}}==

PolynomialInterpolation([[0, 1], [1, 6], [2, 17], [3, 34], [4, 57], [5, 86], [6, 121], [7, 162], [8, 209], [9, 262], [10, 321]], 'x');

Result:

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

 

Using the built-in "Fit" function.

Second version: using built-in "InterpolatingPolynomial" function.

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

 

=={{header|Modula-2}}==

FROM FormatString IMPORT FormatString;

FROM RealStr IMPORT RealToStr;

 

ReadChar;

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

{{output}}

 

 

 

 

 

{{out}}

<pre>

y=3x^2+2x+1

10 321 321

</pre>

 

=={{header|PowerShell}}==

function qr([double[][]]$A) {

$m,$n = $A.count, $A[0].count

"X^2 X constant"

"$(polyfit $x $y 2)"

{{out}}

<pre>

 

{{libheader|NumPy}}

>>> 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|Raku}}==

(formerly Perl 6)

 

 

constant @x1 = <0 1 2 3 4 5 6 7 8 9 10>;

 

constant $y = [+] @y Z* @e;

 

.say for

{{out}}

<pre>1

 

=={{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.t * mx).inv * mx.t * my).transpose.to_a[0].map(&:to_f)

'''Testing:'''

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

{{out}}

<pre>[1.0, 2.0, 3.0]</pre>

{{libheader|Scastie qualified}}

{{works with|Scala|2.13}}

private def xy = Seq(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321).zipWithIndex.map(_.swap)

 

polyRegression(xy)

 

 

=={{header|Sidef}}==

{{trans|Ruby}}

var A = Matrix.build(x.len, degree+1, {|i,j|

x[i]**j

)

 

{{out}}

<pre>[1, 2, 3]</pre>

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

{{trans|Kotlin}}

 

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

 

polyRegression(x: x, y: y)

 

{{out}}

<!-- 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|VBA}}==

Excel VBA has built in capability for line estimation.

Private Function polynomial_regression(y As Variant, x As Variant, degree As Integer) As Variant

Dim a() As Double

Debug.Print "Degrees of freedom:"; result(4, 2)

Debug.Print "Standard error of y estimate:"; result(3, 2)

<pre>coefficients : 1, 2, 3,

standard errors: 0, 0, 0,

Degrees of freedom: 8

Standard error of y estimate: 2,79482284961344E-14 </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>