Polynomial regression: Difference between revisions
→{{header|jq}}
(41 intermediate revisions by 23 users not shown) | |||
Line 11:
This task is intended as a subtask for [[Measure relative performance of sorting algorithms implementations]].
=={{header|11l}}==
{{trans|Swift}}
<syntaxhighlight lang="11l">F average(arr)
R sum(arr) / Float(arr.len)
F poly_regression(x, y)
V xm = average(x)
V ym = average(y)
V x2m = average(x.map(i -> i * i))
V x3m = average(x.map(i -> i ^ 3))
V x4m = average(x.map(i -> i ^ 4))
V xym = average(zip(x, y).map((i, j) -> i * j))
V x2ym = average(zip(x, y).map((i, j) -> i * i * j))
V sxx = x2m - xm * xm
V sxy = xym - xm * ym
V sxx2 = x3m - xm * x2m
V sx2x2 = x4m - x2m * x2m
V sx2y = x2ym - x2m * ym
V b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
V c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
V a = ym - b * xm - c * x2m
F abc(xx)
R (@a + @b * xx) + (@c * xx * xx)
print("y = #. + #.x + #.x^2\n".format(a, b, c))
print(‘ Input Approximation’)
print(‘ x y y1’)
L(i) 0 .< x.len
print(‘#2 #3 #3.1’.format(x[i], y[i], abc(i)))
V x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
V y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]
poly_regression(x, y)</syntaxhighlight>
{{out}}
<pre>
y = 1 + 2x + 3x^2
Input Approximation
x y y1
0 1 1.0
1 6 6.0
2 17 17.0
3 34 34.0
4 57 57.0
5 86 86.0
6 121 121.0
7 162 162.0
8 209 209.0
9 262 262.0
10 321 321.0
</pre>
=={{header|Ada}}==
<
function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is
Line 24 ⟶ 80:
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 ''x''<sub>''i''</sub> and ''y''<sub>''i''</sub>. The basis φ<sub>''j''</sub> is ''x''<sup>''j''</sup>, ''j''=0,1,..,''N''. The implementation is straightforward. First the plane matrix A is created. A<sub>ji</sub>=φ<sub>''j''</sub>(''x''<sub>''i''</sub>). Then the linear problem AA<sup>''T''</sup>''c''=A''y'' is solved. The result ''c''<sub>''j''</sub> 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;
Line 41 ⟶ 97:
Put (C (1), Aft => 3, Exp => 0);
Put (C (2), Aft => 3, Exp => 0);
end Fitting;</
{{out}}
<pre>
Line 54 ⟶ 110:
<!-- {{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
Line 167 ⟶ 223:
);
print polynomial(d)
END # fitting #</
{{out}}
<pre>
3x**2+2x+1
1.0848x**2+10.3552x-0.6164
</pre>
=={{header|AutoHotkey}}==
{{trans|Lua}}
<syntaxhighlight lang="autohotkey">
regression(xa,ya){
n := xa.Count()
xm := ym := x2m := x3m := x4m := xym := x2ym := 0
loop % n {
i := A_Index
xm := xm + xa[i]
ym := ym + ya[i]
x2m := x2m + xa[i] * xa[i]
x3m := x3m + xa[i] * xa[i] * xa[i]
x4m := x4m + xa[i] * xa[i] * xa[i] * xa[i]
xym := xym + xa[i] * ya[i]
x2ym := x2ym + xa[i] * xa[i] * ya[i]
}
xm := xm / n
ym := ym / n
x2m := x2m / n
x3m := x3m / n
x4m := x4m / n
xym := xym / n
x2ym := x2ym / n
sxx := x2m - xm * xm
sxy := xym - xm * ym
sxx2 := x3m - xm * x2m
sx2x2 := x4m - x2m * x2m
sx2y := x2ym - x2m * ym
b := (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
c := (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
a := ym - b * xm - c * x2m
result := "Input`tApproximation`nx y`ty1`n"
loop % n
i := A_Index, result .= xa[i] ", " ya[i] "`t" eval(a, b, c, xa[i]) "`n"
return "y = " c "x^2" " + " b "x + " a "`n`n" result
}
eval(a,b,c,x){
return a + (b + c*x) * x
}</syntaxhighlight>
Examples:<syntaxhighlight lang="autohotkey">xa := [0, 1, 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10]
ya := [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]
MsgBox % result := regression(xa, ya)
return</syntaxhighlight>
{{out}}
<pre>y = 3.000000x^2 + 2.000000x + 1.000000
Input Approximation
x y y1
0, 1 1.000000
1, 6 6.000000
2, 17 17.000000
3, 34 34.000000
4, 57 57.000000
5, 86 86.000000
6, 121 121.000000
7, 162 162.000000
8, 209 209.000000
9, 262 262.000000
10, 321 321.000000</pre>
=={{header|AWK}}==
{{trans|Lua}}
<syntaxhighlight lang="awk">
BEGIN{
i = 0;
xa[i] = 0; i++;
xa[i] = 1; i++;
xa[i] = 2; i++;
xa[i] = 3; i++;
xa[i] = 4; i++;
xa[i] = 5; i++;
xa[i] = 6; i++;
xa[i] = 7; i++;
xa[i] = 8; i++;
xa[i] = 9; i++;
xa[i] = 10; i++;
i = 0;
ya[i] = 1; i++;
ya[i] = 6; i++;
ya[i] = 17; i++;
ya[i] = 34; i++;
ya[i] = 57; i++;
ya[i] = 86; i++;
ya[i] =121; i++;
ya[i] =162; i++;
ya[i] =209; i++;
ya[i] =262; i++;
ya[i] =321; i++;
exit;
}
{
# (nothing to do)
}
END{
a = 0; b = 0; c = 0; # globals - will change by regression()
regression(xa,ya);
printf("y = %6.2f x^2 + %6.2f x + %6.2f\n",c,b,a);
printf("%-13s %-8s\n","Input","Approximation");
printf("%-6s %-6s %-8s\n","x","y","y^")
for (i=0;i<length(xa);i++) {
printf("%6.1f %6.1f %8.3f\n",xa[i],ya[i],eval(a,b,c,xa[i]));
}
}
function eval(a,b,c,x) {
return a+b*x+c*x*x;
}
# locals
function regression(x,y, n,xm,ym,x2m,x3m,x4m,xym,x2ym,sxx,sxy,sxx2,sx2x2,sx2y) {
n = 0
xm = 0.0;
ym = 0.0;
x2m = 0.0;
x3m = 0.0;
x4m = 0.0;
xym = 0.0;
x2ym = 0.0;
for (i in x) {
xm += x[i];
ym += y[i];
x2m += x[i] * x[i];
x3m += x[i] * x[i] * x[i];
x4m += x[i] * x[i] * x[i] * x[i];
xym += x[i] * y[i];
x2ym += x[i] * x[i] * y[i];
n++;
}
xm = xm / n;
ym = ym / n;
x2m = x2m / n;
x3m = x3m / n;
x4m = x4m / n;
xym = xym / n;
x2ym = x2ym / n;
sxx = x2m - xm * xm;
sxy = xym - xm * ym;
sxx2 = x3m - xm * x2m;
sx2x2 = x4m - x2m * x2m;
sx2y = x2ym - x2m * ym;
b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
a = ym - b * xm - c * x2m;
}
</syntaxhighlight>
{{out}}
<pre>
y = 3.00 x^2 + 2.00 x + 1.00
Input Approximation
x y y^
0.0 1.0 1.000
1.0 6.0 6.000
2.0 17.0 17.000
3.0 34.0 34.000
4.0 57.0 57.000
5.0 86.0 86.000
6.0 121.0 121.000
7.0 162.0 162.000
8.0 209.0 209.000
9.0 262.0 262.000
10.0 321.0 321.000
</pre>
Line 179 ⟶ 406:
and fits an order-5 polynomial, so the test data for this task
is hardly challenging!
<
Max% = 10000
Line 227 ⟶ 454:
FOR term% = 5 TO 0 STEP -1
PRINT ;vector(term%) " * x^" STR$(term%)
NEXT</
{{out}}
<pre>
Line 243 ⟶ 470:
'''Include''' file (to make the code reusable easily) named <tt>polifitgsl.h</tt>
<
#define _POLIFITGSL_H
#include <gsl/gsl_multifit.h>
Line 250 ⟶ 477:
bool polynomialfit(int obs, int degree,
double *dx, double *dy, double *store); /* n, p */
#endif</
'''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,
Line 292 ⟶ 519:
return true; /* we do not "analyse" the result (cov matrix mainly)
to know if the fit is "good" */
}</
'''Testing''':
<
#include "polifitgsl.h"
Line 314 ⟶ 541:
}
return 0;
}</
{{out}}
<pre>1.000000
Line 322 ⟶ 549:
=={{header|C sharp|C#}}==
{{libheader|Math.Net}}
<
{
// Vandermonde matrix
Line 336 ⟶ 563:
var p = r.Inverse().Multiply(q.TransposeThisAndMultiply(yv));
return p.Column(0).ToArray();
}</
Example:
<
{
const int degree = 2;
Line 349 ⟶ 576:
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>
Line 414 ⟶ 641:
return 0;
}</
{{out}}
<pre>y = 1 + 2x + 3x^2
Line 434 ⟶ 661:
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)))
Line 442 ⟶ 669:
(setf (aref A i j)
(expt (aref x 0 i) j))))
(lsqr A (mtp y))))</
Example:
<
(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))</
=={{header|D}}==
{{trans|Kotlin}}
<
import std.range;
import std.stdio;
Line 500 ⟶ 727:
auto y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
polyRegression(x, y);
}</
{{out}}
<pre>y = 1 + 2x + 3x^2
Line 516 ⟶ 743:
9 262 262.0
10 321 321.0</pre>
=={{header|EasyLang}}==
{{trans|Lua}}
<syntaxhighlight lang=easylang>
func eval a b c x .
return a + (b + c * x) * x
.
proc regression xa[] ya[] . .
n = len xa[]
for i = 1 to n
xm = xm + xa[i]
ym = ym + ya[i]
x2m = x2m + xa[i] * xa[i]
x3m = x3m + xa[i] * xa[i] * xa[i]
x4m = x4m + xa[i] * xa[i] * xa[i] * xa[i]
xym = xym + xa[i] * ya[i]
x2ym = x2ym + xa[i] * xa[i] * ya[i]
.
xm = xm / n
ym = ym / n
x2m = x2m / n
x3m = x3m / n
x4m = x4m / n
xym = xym / n
x2ym = x2ym / n
#
sxx = x2m - xm * xm
sxy = xym - xm * ym
sxx2 = x3m - xm * x2m
sx2x2 = x4m - x2m * x2m
sx2y = x2ym - x2m * ym
#
b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
a = ym - b * xm - c * x2m
print "y = " & a & " + " & b & "x + " & c & "x^2"
numfmt 0 3
for i = 1 to n
print xa[i] & " " & ya[i] & " " & eval a b c xa[i]
.
.
xa[] = [ 0 1 2 3 4 5 6 7 8 9 10 ]
ya[] = [ 1 6 17 34 57 86 121 162 209 262 321 ]
regression xa[] ya[]
</syntaxhighlight>
=={{header|Emacs Lisp}}==
{{libheader|Calc}}
<syntaxhighlight lang="lisp">(let ((x '(0 1 2 3 4 5 6 7 8 9 10))
(y '(1 6 17 34 57 86 121 162 209 262 321)))
(calc-eval "fit(a*x^2+b*x+c,[x],[a,b,c],[$1 $2])" nil (cons 'vec x) (cons 'vec y)))</syntaxhighlight>
{{out}}
"3. x^2 + 1.99999999996 x + 1.00000000006"
=={{header|Fortran}}==
{{libheader|LAPACK}}
<
contains
Line 598 ⟶ 865:
end function
end module</
===Example===
<
use fitting
implicit none
Line 619 ⟶ 886:
write (*, '(F9.4)') a
end program</
{{out}} (lower powers first, so this seems the opposite of the Python output):
Line 629 ⟶ 896:
=={{header|FreeBASIC}}==
General regressions for different polynomials, here it is for degree 2, (3 terms).
<syntaxhighlight lang="freebasic">#Include "crt.bi" 'for rounding only
Type vector
Dim As Double element(Any)
End Type
Type
Dim As Double
Declare Function inverse() As matrix
Declare Function transpose() As matrix
private:
Declare Function GaussJordan(As vector) As vector
End Type
'mult operators
Operator *(m1 As matrix,m2 As matrix) As matrix
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,2)<>Ubound(m2.element,1) Then
Print "Can't do"
Exit Operator
End If
Dim As matrix ans
Redim ans.element(rows,columns)
Dim rxc As Double
For r As Integer=1 To rows
For c As Integer=1 To columns
rxc=0
For k As Integer = 1 To Ubound(m1.element,2)
rxc=rxc+m1.element(r,k)*m2.element(k,c)
Next k
ans.element(r,c)=rxc
Next c
Next r
Operator= ans
End Operator
Operator *(m1 As matrix,m2 As vector) As vector
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,2)<>Ubound(m2.element) Then
Print "Can't do"
Exit Operator
End If
Dim As vector ans
Redim ans.element(rows)
Dim rxc As Double
For r As Integer=1 To rows
rxc=0
For k As Integer = 1 To Ubound(m1.element,2)
rxc=rxc+m1.element(r,k)*m2.element(k)
Next k
ans.element(r)=rxc
Next r
Operator= ans
End Operator
Function matrix.transpose() As matrix
Dim As matrix b
Redim b.element(1 To Ubound(this.element,2),1 To Ubound(this.element,1))
For i As Long=1 To Ubound(this.element,1)
For j As Long=1 To Ubound(this.element,2)
b.element(j,i)=this.element(i,j)
Next
Next
Return b
End Function
Function matrix.GaussJordan(rhs As vector) As vector
Dim As Integer n=Ubound(rhs.element)
Dim As vector ans=rhs,r=rhs
Dim As matrix b=This
#macro pivot(num)
For p1 As Integer = num To n - 1
For p2 As Integer = p1 + 1 To n
If Abs(b.element(p1,num))<Abs(b.element(p2,num)) Then
Swap r.element(p1),r.element(p2)
For g As Integer=1 To n
Swap b.element(p1,g),b.element(p2,g)
Next g
End If
Next p2
Next p1
#endmacro
For k As Integer=1 To n-1
pivot(k)
For row As Integer =k To n-1
If b.element(row+1,k)=0 Then Exit For
Var f=b.element(k,k)/b.element(row+1,k)
r.element(row+1)=r.element(row+1)*f-r.element(k)
For g As Integer=1 To n
b.element((row+1),g)=b.element((row+1),g)*f-b.element(k,g)
Next g
Next row
Next k
'back substitute
For z As Integer=n To 1 Step -1
ans.element(z)=r.element(z)/b.element(z,z)
For j As Integer = n To z+1 Step -1
ans.element(z)=ans.element(z)-(b.element(z,j)*ans.element(j)/b.element(z,z))
Next j
End Function
Function matrix.inverse() As matrix
Var ub1=Ubound(this.element,1),ub2=Ubound(this.element,2)
Dim As matrix
Dim As vector
Redim temp.element(1 To ub1):Redim null_.element(1 To ub1)
Redim ans.element(1 To ub1,1 To ub2)
For a As
ans.element(b,a)=temp.element(b)
Next b
Next a
Return ans
End Function
'vandermode of x
Function vandermonde(x_values() As Double,w As Long) As matrix
Dim As matrix mat
Var n=Ubound(x_values)
Redim mat.element(1 To n,1 To w)
For a As Integer=1 To n
For b As Integer=1 To w
mat.element(a,b)=x_values(a)^(b-1)
Next b
Next a
Return mat
End Function
'main preocedure
Sub regress(x_values() As Double,y_values() As Double,ans() As Double,n As Long)
Redim ans(1 To Ubound(x_values))
Dim As matrix m1= vandermonde(x_values(),n)
Dim As matrix T=m1.transpose
Dim As vector y
Redim y.element(1 To Ubound(ans))
For n As Long=1 To Ubound(y_values)
y.element(n)=y_values(n)
Next n
Dim As vector result=(((T*m1).inverse)*T)*y
Redim Preserve ans(1 To n)
For n As Long=1 To Ubound(ans)
ans(n)=result.element(n)
Next n
End Sub
'Evaluate a polynomial at x
Function polyeval(Coefficients() As Double,Byval x As Double) As Double
Dim As Double acc
For i As Long=Ubound(Coefficients) To Lbound(Coefficients) Step -1
acc=acc*x+Coefficients(i)
Next i
Return acc
End Function
Function CRound(Byval x As Double,Byval precision As Integer=30) As String
If precision>30 Then precision=30
Dim As zstring * 40 z:Var s="%." &str(Abs(precision)) &"f"
sprintf(z,s,x)
If Val(z) Then Return Rtrim(Rtrim(z,"0"),".")Else Return "0"
End Function
Function show(a() As Double,places as long=10) As String
Dim As String s,g
For n As Long=Lbound(a) To Ubound(a)
If n<3 Then g="" Else g="^"+Str(n-1)
if val(cround(a(n),places))<>0 then
s+= Iif(Sgn(a(n))>=0,"+","")+cround(a(n),places)+ Iif(n=Lbound(a),"","*x"+g)+" "
end if
Next n
Return s
End Function
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 ans()
regress(x(),y(),ans(),3)
print show(ans())
sleep</syntaxhighlight>
{{out}}
<pre>+1 +2*x +3*x^2</pre>
=={{header|GAP}}==
<
local a;
a := List([0 .. n], i -> List(x, s -> s^i));
Line 734 ⟶ 1,100:
# Return coefficients in ascending degree order
PolynomialRegression(x, y, 2);
# [ 1, 2, 3 ]</
=={{header|gnuplot}}==
<
f(x) = a*x**2 + b*x + c
Line 761 ⟶ 1,127:
e
print sprintf("\n --- \n Polynomial fit: %.4f x^2 + %.4f x + %.4f\n", a, b, c)</
=={{header|Go}}==
===Library gonum/matrix===
<
import (
"log"
"gonum.org/v1/gonum/mat"
)
func main() {
var (
y = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
a = Vandermonde(x, degree+1)
b = mat.NewDense(len(y), 1, y)
c
)
const trans = false
err := qr.SolveTo(c, trans, b)
if err != nil {
log.Fatalf("could not solve QR: %+v", err)
}
fmt.Printf("%.3f\n", mat.Formatted(c))
}
func Vandermonde(a []float64,
}
}
}</
{{out}}
<pre>
Line 814 ⟶ 1,181:
===Library go.matrix===
Least squares solution using QR decomposition and package [http://github.com/skelterjohn/go.matrix go.matrix].
<
import (
Line 854 ⟶ 1,221:
}
fmt.Println(c)
}</
{{out}} (lowest order coefficient first)
<pre>
Line 862 ⟶ 1,229:
=={{header|Haskell}}==
Uses module Matrix.LU from [http://hackage.haskell.org/package/dsp hackageDB DSP]
<
import Data.Array
import Control.Monad
Line 872 ⟶ 1,239:
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</
{{out}} in GHCi:
<
[1.0,2.0,3.0]</
=={{header|HicEst}}==
<
x = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
Line 892 ⟶ 1,259:
! 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</
{{out}}
<pre>SOLVE performs a (nonlinear) least-square fit (Levenberg-Marquardt):
Line 898 ⟶ 1,265:
=={{header|Hy}}==
<
(setv x (range 11))
(setv y [1 6 17 34 57 86 121 162 209 262 321])
(print (polyfit x y 2))</
=={{header|J}}==
<
(%. ^/~@x:@i.@#) Y
1 2 3 0 0 0 0 0 0 0 0</
Note that this implementation does not use floating point numbers,
Line 916 ⟶ 1,283:
but for small problems like this it is inconsequential.
The above solution fits a polynomial of order 11 (or, more specifically, a polynomial whose order matches the length of its argument sequence).
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:
<
1 2 3</
(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.)
Line 927 ⟶ 1,294:
{{trans|D}}
{{works with|Java|8}}
<
import java.util.function.IntToDoubleFunction;
import java.util.stream.IntStream;
Line 934 ⟶ 1,301:
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 x2m = Arrays.stream(
double x3m = Arrays.stream(
double x4m = Arrays.stream(
double xym = 0.0;
for (int i = 0; i < x.length && i < y.length; ++i) {
Line 976 ⟶ 1,342:
polyRegression(x, y);
}
}</
{{out}}
<pre>y = 1.0 + 2.0x + 3.0x^2
Line 992 ⟶ 1,358:
9 262 262.0
10 321 321.0</pre>
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
'''Works with jq, the C implementation of jq'''
'''Works with gojq, the Go implementation of jq'''
'''Works with jaq, the Rust implementation of jq'''
<syntaxhighlight lang="jq">
def mean: add/length;
def inner_product($y):
. as $x
| reduce range(0; length) as $i (0; . + ($x[$i] * $y[$i]));
# $x and $y should be arrays of the same length
# Emit { a, b, c, z}
# Attempt to avoid overflow
def polynomialRegression($x; $y):
($x | length) as $length
| ($length * $length) as $l2
| ($x | map(./$length)) as $xs
| ($xs | add) as $xm
| ($y | mean) as $ym
| ($xs | map(. * .) | add * $length) as $x2m
| ($x | map( (./$length) * . * .) | add) as $x3m
| ($xs | map(. * . | (.*.) ) | add * $l2 * $length) as $x4m
| ($xs | inner_product($y)) as $xym
| ($xs | map(. * .) | inner_product($y) * $length) as $x2ym
| ($x2m - $xm * $xm) as $sxx
| ($xym - $xm * $ym) as $sxy
| ($x3m - $xm * $x2m) as $sxx2
| ($x4m - $x2m * $x2m) as $sx2x2
| ($x2ym - $x2m * $ym) as $sx2y
| {z: ([$x,$y] | transpose) }
| .b = ($sxy * $sx2x2 - $sx2y * $sxx2) / ($sxx * $sx2x2 - $sxx2 * $sxx2)
| .c = ($sx2y * $sxx - $sxy * $sxx2) / ($sxx * $sx2x2 - $sxx2 * $sxx2)
| .a = $ym - .b * $xm - .c * $x2m ;
# Input: {a,b,c,z}
def report:
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
def abc($x): .a + .b * $x + .c * $x * $x;
def print($p): "\($p[0] | lpad(3)) \($p[1] | lpad(4)) \(abc($p[0]) | lpad(5))";
"y = \(.a) + \(.b)x + \(.c)x^2\n",
" Input Approximation",
" x y y\u0302",
print(.z[]) ;
def x: [range(0;11)];
def y: [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
polynomialRegression(x; y)
| report
</syntaxhighlight>
{{output}}
<pre>
y = 1 + 2x + 3x^2
Input Approximation
x y ŷ
0 1 1
1 6 6
2 17 17
3 34 34
4 57 57
5 86 86
6 121 121
7 162 162
8 209 209
9 262 262
10 321 321
</pre>
=={{header|Julia}}==
Line 997 ⟶ 1,439:
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]
@show polyfit(x, y, 2)</
{{out}}
Line 1,008 ⟶ 1,450:
=={{header|Kotlin}}==
{{trans|REXX}}
<
fun polyRegression(x: IntArray, y: IntArray) {
val xm = x.average()
val ym = y.average()
val x2m =
val x3m =
val x4m =
val xym = x.zip(y).map { it.first * it.second }.average()
val x2ym = x.zip(y).map { it.first * it.first * it.second }.average()
Line 1,036 ⟶ 1,476:
println(" Input Approximation")
println(" x y y1")
for (
System.out.printf("%2d %3d %5.1f\n",
}
}
fun main(
val x = IntArray(11) { it }
val y = intArrayOf(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
polyRegression(x, y)
}</
{{out}}
Line 1,068 ⟶ 1,508:
=={{header|Lua}}==
{{trans|Modula-2}}
<
return a + (b + c * x) * x
end
Line 1,119 ⟶ 1,559:
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}
regression(xa, ya)</
{{out}}
<pre>y = 1 + 2x + 3x^2
Line 1,135 ⟶ 1,575:
=={{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');
</syntaxhighlight>
Result:
<pre>3*x^2+2*x+1</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
Using the built-in "Fit" function.
<syntaxhighlight 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]</syntaxhighlight>
Second version: using built-in "InterpolatingPolynomial" function.
<
WolframAlpha version:
<syntaxhighlight lang="mathematica">curve fit (0,1), (1,6), (2,17), (3,34), (4,57), (5,86), (6,121), (7,162), (8,209), (9,262), (10,321)</syntaxhighlight>
Result:
<pre>1 + 2x + 3x^2</pre>
Line 1,157 ⟶ 1,597:
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)
Line 1,163 ⟶ 1,603:
ans =
2.999999999999998 2.000000000000019 0.999999999999956</
=={{header|МК-61/52}}==
Part 1:
<syntaxhighlight lang="text">П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'': В/О x<sub>1</sub> С/П y<sub>1</sub> С/П x<sub>2</sub> С/П y<sub>2</sub> С/П ...
Part 2:
<syntaxhighlight lang="text">ИП5 ПC ИП6 ПD П2 ИП7 П3 ИП4 ИПD *
ИПC ИП5 * - ПD ИП4 ИП7 * ИПC ИП6
* - П7 ИП4 ИПA * ИПC ИП9 * -
Line 1,184 ⟶ 1,624:
ИПD ИП8 * ИП7 ИП3 * - / ПB ИПA
ИПB ИП7 * - ИПD / ПA ИП9 ИПB ИП6
* - ИПA ИП5 * - ИП4 / П9 С/П</
''Result'': Р9 = a<sub>0</sub>, РA = a<sub>1</sub>, РB = a<sub>2</sub>.
=={{header|Modula-2}}==
<
FROM FormatString IMPORT FormatString;
FROM RealStr IMPORT RealToStr;
Line 1,275 ⟶ 1,715:
ReadChar;
END PolynomialRegression.</
=={{header|
{{trans|Kotlin}}
<syntaxhighlight lang="nim">import lenientops, sequtils, stats, strformat
proc polyRegression(x, y: openArray[int]) =
let xm = mean(x)
let ym = mean(y)
let x2m = mean(x.mapIt(it * it))
let x3m = mean(x.mapIt(it * it * it))
let x4m = mean(x.mapIt(it * it * it * it))
let xym = mean(zip(x, y).mapIt(it[0] * it[1]))
let x2ym = mean(zip(x, y).mapIt(it[0] * it[0] * it[1]))
let sxx = x2m - xm * xm
let sxy = xym - xm * ym
let sxx2 = x3m - xm * x2m
let sx2x2 = x4m - x2m * x2m
let sx2y = x2ym - x2m * ym
let b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
let c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
let a = ym - b * xm - c * x2m
func abc(x: int): float = a + b * x + c * x * x
echo &"y = {a} + {b}x + {c}x²\n"
echo " Input Approximation"
echo " x y y1"
for (xi, yi) in zip(x, y):
echo &"{xi:2} {yi:3} {abc(xi):5}"
let x = toSeq(0..10)
let y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]
polyRegression(x, y)</syntaxhighlight>
{{out}}
<pre>y = 1.0 + 2.0x + 3.0x²
Input Approximation
x y y1
0 1 1
1 6 6
2 17 17
3 34 34
4 57 57
5 86 86
6 121 121
7 162 162
8 209 209
9 262 262
10 321 321</pre>
=={{header|OCaml}}==
{{trans|Kotlin}}
{{libheader|Base}}
<syntaxhighlight lang="ocaml">open Base
open Stdio
let mean fa =
let open Float in
(Array.reduce_exn fa ~f:(+)) / (of_int (Array.length fa))
let regression xs ys =
let open Float in
let xm = mean xs in
let ym = mean ys in
let x2m = Array.map xs ~f:(fun x -> x * x) |> mean in
let x3m = Array.map xs ~f:(fun x -> x * x * x) |> mean in
let x4m = Array.map xs ~f:(fun x -> let x2 = x * x in x2 * x2) |> mean in
let xzipy = Array.zip_exn xs ys in
let xym = Array.map xzipy ~f:(fun (x, y) -> x * y) |> mean in
let x2ym = Array.map xzipy ~f:(fun (x, y) -> x * x * y) |> mean in
let sxx = x2m - xm * xm in
let sxy = xym - xm * ym in
let sxx2 = x3m - xm * x2m in
let sx2x2 = x4m - x2m * x2m in
let sx2y = x2ym - x2m * ym in
let b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) in
let c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) in
let a = ym - b * xm - c * x2m in
let abc xx = a + b * xx + c * xx * xx in
printf "y = %.1f + %.1fx + %.1fx^2\n\n" a b c;
printf " Input Approximation\n";
printf " x y y1\n";
Array.iter xzipy ~f:(fun (xi, yi) ->
printf "%2g %3g %5.1f\n" xi yi (abc xi)
)
let () =
let x = Array.init 11 ~f:Float.of_int in
let y = [| 1.; 6.; 17.; 34.; 57.; 86.; 121.; 162.; 209.; 262.; 321. |] in
regression x y</syntaxhighlight>
{{out}}
<pre>
y = 1.0 + 2.0x + 3.0x^2
Input Approximation
x y y1
0 1 1.0
1 6 6.0
2 17 17.0
3 34 34.0
4 57 57.0
5 86 86.0
6 121 121.0
7 162 162.0
8 209 209.0
9 262 262.0
10 321 321.0
</pre>
=={{header|Octave}}==
<syntaxhighlight lang="octave">x = [0:10];
y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
coeffs = polyfit(x, y, 2)</syntaxhighlight>
=={{header|PARI/GP}}==
Lagrange interpolating polynomial:
<syntaxhighlight 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])</syntaxhighlight>
In newer versions, this can be abbreviated:
<syntaxhighlight lang="parigp">polinterpolate([0..10],[1,6,17,34,57,86,121,162,209,262,321])</syntaxhighlight>
{{out}}
<pre>3*x^2 + 2*x + 1</pre>
Least-squares fit:
<syntaxhighlight 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))</syntaxhighlight>
<small>Code thanks to [http://pari.math.u-bordeaux.fr/archives/pari-users-1105/msg00006.html Bill Allombert]</small>
{{out}}
<pre>3*x^2 + 2*x + 1</pre>
Least-squares polynomial fit in its own function:
<syntaxhighlight lang="parigp">lsf(X,Y,n)=my(M=matrix(#X,n+1,i,j,X[i]^(j-1))); Polrev(matsolve(M~*M,M~*Y~))
lsf([0..10], [1,6,17,34,57,86,121,162,209,262,321], 2)</syntaxhighlight>
=={{header|Perl}}==
This code identical to that of [[Multiple regression]] task.
<syntaxhighlight lang="perl">use strict;
use warnings;
use Statistics::Regression;
my @x = <0 1 2 3 4 5 6 7 8 9 10>;
my @y = <1 6 17 34 57 86 121 162 209 262 321>;
my @model = ('const', 'X', 'X**2');
my $reg = Statistics::Regression->new( '', [@model] );
$reg->include( $y[$_], [ 1.0, $x[$_], $x[$_]**2 ]) for 0..@y-1;
my @coeff = $reg->theta();
printf "%-6s %8.3f\n", $model[$_], $coeff[$_] for 0..@model-1;</syntaxhighlight>
{{output}}
<pre>const 1.000
X**2 3.000</pre>
PDL Alternative:
<syntaxhighlight lang="perl">#!/usr/bin/perl -w
use strict;
use PDL;
use PDL::Math;
use PDL::Fit::Polynomial;
my $x = float [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
my $y = float [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
# above will output: 3.00000037788248 * $x**2 + 1.99999750988868 * $x + 1.00000180493936
# $x = float [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
# $y = float [ 2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0];
# above correctly returns: " 1.08484845125187 * $x**2 + 10.3551513321297 * $x-0.616363852007752 "
my ($yfit, $coeffs) = fitpoly1d $x, $y, 3; # 3rd degree
foreach (reverse(0..$coeffs->dim(0)-1)) {
print " +" unless(($coeffs->at($_) <0) || $_==$coeffs->dim(0)-1); # let the unary minus replace the + operator
print " ";
print $coeffs->at($_);
print " * \$x" if($_);
print "**$_" if($_>1);
print "\n" unless($_)
}
</syntaxhighlight>
{{output}}
<pre> 3.00000037788248 * $x**2 + 1.99999750988868 * $x + 1.00000180493936</pre>
=={{header|Phix}}==
{{trans|REXX}}
{{libheader|Phix/online}}
{{libheader|Phix/pGUI}}
You can run this online [http://phix.x10.mx/p2js/Polynomial_regression.htm here].
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Polynomial_regression.exw</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span><span style="color: #000000;">7</span><span style="color: #0000FF;">,</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span><span style="color: #000000;">10</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">6</span><span style="color: #0000FF;">,</span><span style="color: #000000;">17</span><span style="color: #0000FF;">,</span><span style="color: #000000;">34</span><span style="color: #0000FF;">,</span><span style="color: #000000;">57</span><span style="color: #0000FF;">,</span><span style="color: #000000;">86</span><span style="color: #0000FF;">,</span><span style="color: #000000;">121</span><span style="color: #0000FF;">,</span><span style="color: #000000;">162</span><span style="color: #0000FF;">,</span><span style="color: #000000;">209</span><span style="color: #0000FF;">,</span><span style="color: #000000;">262</span><span style="color: #0000FF;">,</span><span style="color: #000000;">321</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">regression</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">xm</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ym</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x2m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x3m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x4m</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">xym</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">x2ym</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">xi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">yi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">y</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">xm</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">xi</span>
<span style="color: #000000;">ym</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">yi</span>
<span style="color: #000000;">x2m</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xi</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">x3m</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xi</span><span style="color: #0000FF;">,</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">x4m</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xi</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">xym</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">xi</span><span style="color: #0000FF;">*</span><span style="color: #000000;">yi</span>
<span style="color: #000000;">x2ym</span> <span style="color: #0000FF;">+=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xi</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">yi</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">xm</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">n</span>
<span style="color: #000000;">ym</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">n</span>
<span style="color: #000000;">x2m</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">n</span>
<span style="color: #000000;">x3m</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">n</span>
<span style="color: #000000;">x4m</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">n</span>
<span style="color: #000000;">xym</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">n</span>
<span style="color: #000000;">x2ym</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">n</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">Sxx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x2m</span><span style="color: #0000FF;">-</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">xm</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">Sxy</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">xym</span><span style="color: #0000FF;">-</span><span style="color: #000000;">xm</span><span style="color: #0000FF;">*</span><span style="color: #000000;">ym</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">Sxx2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x3m</span><span style="color: #0000FF;">-</span><span style="color: #000000;">xm</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x2m</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">Sx2x2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x4m</span><span style="color: #0000FF;">-</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x2m</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">Sx2y</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">x2ym</span><span style="color: #0000FF;">-</span><span style="color: #000000;">x2m</span><span style="color: #0000FF;">*</span><span style="color: #000000;">ym</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">B</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">Sxy</span><span style="color: #0000FF;">*</span><span style="color: #000000;">Sx2x2</span><span style="color: #0000FF;">-</span><span style="color: #000000;">Sx2y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">Sxx2</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">Sxx</span><span style="color: #0000FF;">*</span><span style="color: #000000;">Sx2x2</span><span style="color: #0000FF;">-</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">Sxx2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">C</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">Sx2y</span><span style="color: #0000FF;">*</span><span style="color: #000000;">Sxx</span><span style="color: #0000FF;">-</span><span style="color: #000000;">Sxy</span><span style="color: #0000FF;">*</span><span style="color: #000000;">Sxx2</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">Sxx</span><span style="color: #0000FF;">*</span><span style="color: #000000;">Sx2x2</span><span style="color: #0000FF;">-</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">Sxx2</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">A</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ym</span><span style="color: #0000FF;">-</span><span style="color: #000000;">B</span><span style="color: #0000FF;">*</span><span style="color: #000000;">xm</span><span style="color: #0000FF;">-</span><span style="color: #000000;">C</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x2m</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">C</span><span style="color: #0000FF;">,</span><span style="color: #000000;">B</span><span style="color: #0000FF;">,</span><span style="color: #000000;">A</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">regression</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">b</span><span style="color: #0000FF;">*</span><span style="color: #000000;">x</span><span style="color: #0000FF;">+</span><span style="color: #000000;">c</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"y=%gx^2+%gx+%g\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span><span style="color: #000000;">c</span><span style="color: #0000FF;">})</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n x y f(x)\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" %2d %3d %3g\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">y</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000080;font-style:italic;">-- And a simple plot (re-using x,y from above)</span>
<span style="color: #008080;">include</span> <span style="color: #000000;">pGUI</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #008080;">include</span> <span style="color: #000000;">IupGraph</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">get_data</span><span style="color: #0000FF;">(</span><span style="color: #004080;">Ihandle</span> <span style="color: #000000;">graph</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">w</span><span style="color: #0000FF;">,</span><span style="color: #000000;">h</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">IupGetIntInt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">graph</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"DRAWSIZE"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">IupSetInt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">graph</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"YTICK"</span><span style="color: #0000FF;">,</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;"><</span><span style="color: #000000;">240</span><span style="color: #0000FF;">?</span><span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">h</span><span style="color: #0000FF;"><</span><span style="color: #000000;">150</span><span style="color: #0000FF;">?</span><span style="color: #000000;">80</span><span style="color: #0000FF;">:</span><span style="color: #000000;">40</span><span style="color: #0000FF;">):</span><span style="color: #000000;">20</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">,</span><span style="color: #004600;">CD_RED</span><span style="color: #0000FF;">}}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #7060A8;">IupOpen</span><span style="color: #0000FF;">()</span>
<span style="color: #004080;">Ihandle</span> <span style="color: #000000;">graph</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">IupGraph</span><span style="color: #0000FF;">(</span><span style="color: #000000;">get_data</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"RASTERSIZE=640x440"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">IupSetAttributes</span><span style="color: #0000FF;">(</span><span style="color: #000000;">graph</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"XTICK=1,XMIN=0,XMAX=10"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">IupSetAttributes</span><span style="color: #0000FF;">(</span><span style="color: #000000;">graph</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"YTICK=20,YMIN=0,YMAX=320"</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">Ihandle</span> <span style="color: #000000;">dlg</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">IupDialog</span><span style="color: #0000FF;">(</span><span style="color: #000000;">graph</span><span style="color: #0000FF;">,</span><span style="color: #008000;">`TITLE="simple plot"`</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">IupSetAttributes</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dlg</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"MINSIZE=245x150"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">IupShow</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dlg</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()!=</span><span style="color: #004600;">JS</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">IupMainLoop</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">IupClose</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<!--</syntaxhighlight>-->
{{out}}
(plus a simple graphical plot, as per [[Polynomial_regression#Racket|Racket]])
<pre>
y=3x^2+2x+1
Line 1,501 ⟶ 2,007:
10 321 321
</pre>
=={{header|PowerShell}}==
<syntaxhighlight lang="powershell">
function qr([double[][]]$A) {
$m,$n = $A.count, $A[0].count
Line 1,609 ⟶ 2,090:
"X^2 X constant"
"$(polyfit $x $y 2)"
</syntaxhighlight>
{{out}}
<pre>
Line 1,620 ⟶ 2,101:
{{libheader|NumPy}}
<
>>> 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
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.])</
Find max absolute error:
<
'1e-013'</
===Example===
For input arrays `x' and `y':
<
>>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0]</
<
>>> 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
Line 1,653 ⟶ 2,134:
which will find the least squares solution via a QR decomposition:
<syntaxhighlight 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)))</
{{out}}
Line 1,666 ⟶ 2,147:
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)
Line 1,690 ⟶ 2,171:
(plot (list (points (map vector xs ys))
(function (poly (fit xs ys 2)))))
</syntaxhighlight>
{{out}}
[[File:polyreg-racket.png]]
Line 1,696 ⟶ 2,177:
=={{header|Raku}}==
(formerly Perl 6)
We'll use a Clifford algebra library. Very slow.
Rationale (in French for some reason):
Le système d'équations peut s'écrire :
<math>\left(a + b x_i + cx_i^2 = y_i\right)_{i=1\ldots N}</math>, où on cherche <math>(a,b,c)\in\mathbb{R}^3</math>. On considère <math>\mathbb{R}^N</math> et on répartit chaque équation sur chaque dimension:
<math> (a + b x_i + cx_i^2)\mathbf{e}_i = y_i\mathbf{e}_i</math>
Posons alors :
<math>
\mathbf{x}_0 = \sum_{i=1}^N \mathbf{e}_i,\,
\mathbf{x}_1 = \sum_{i=1}^N x_i\mathbf{e}_i,\,
\mathbf{x}_2 = \sum_{i=1}^N x_i^2\mathbf{e}_i,\,
\mathbf{y} = \sum_{i=1}^N y_i\mathbf{e}_i
</math>
Le système d'équations devient : <math>a\mathbf{x}_0+b\mathbf{x}_1+c\mathbf{x}_2 = \mathbf{y}</math>.
D'où :
<math>\begin{align}
a = \mathbf{y}\and\mathbf{x}_1\and\mathbf{x}_2/(\mathbf{x}_0\and\mathbf{x_1}\and\mathbf{x_2})\\
b = \mathbf{y}\and\mathbf{x}_2\and\mathbf{x}_0/(\mathbf{x}_1\and\mathbf{x_2}\and\mathbf{x_0})\\
c = \mathbf{y}\and\mathbf{x}_0\and\mathbf{x}_1/(\mathbf{x}_2\and\mathbf{x_0}\and\mathbf{x_1})\\
\end{align}</math>
<syntaxhighlight lang="raku" line>use MultiVector;
constant @x1 = <0 1 2 3 4 5 6 7 8 9 10>;
Line 1,708 ⟶ 2,214:
constant $y = [+] @y Z* @e;
.say for
$y∧$x1∧$x2/($x0∧$x1∧$x2),
$y∧$x2∧$x0/($x1∧$x2∧$x0),
$y∧$x0∧$x1/($x2∧$x0∧$x1);
</syntaxhighlight>
{{out}}
<pre>1
Line 1,725 ⟶ 2,227:
=={{header|REXX}}==
<
* Implementation of http://keisan.casio.com/exec/system/14059932254941
*--------------------------------------------------------------------*/
Line 1,775 ⟶ 2,277:
fun:
Parse Arg x
Return a+b*x+c*x**2 </
{{out}}
<pre>y=1+2*x+3*x**2
Line 1,791 ⟶ 2,293:
9 262 262.000
10 321 321.000</pre>
=={{header|RPL}}==
{{trans|Ada}}
≪ 1 + → x y n
≪ { } n + x SIZE + 0 CON
1 x SIZE '''FOR''' j
1 n '''FOR''' k
{ } k + j + x j GET k 1 - ^ PUT
'''NEXT NEXT'''
DUP y * SWAP DUP TRN * /
<span style="color:grey">@ the following lines convert the resulting vector into a polynomial equation</span>
DUP 'x' STO 1 GET
2 x SIZE '''FOR''' j 'X' * x j GET + '''NEXT'''
EXPAN COLCT
≫ ≫ '<span style="color:blue">FIT</span>' STO
[1 2 3 4 5 6 7 8 9 10] [1 6 17 34 57 86 121 162 209 262 321] 2 <span style="color:blue">FIT</span>
{{out}}
<pre>
1: '3+2*X+1*X^2'
</pre>
=={{header|Ruby}}==
<
def regress x, y, degree
Line 1,802 ⟶ 2,325:
((mx.t * mx).inv * mx.t * my).transpose.to_a[0].map(&:to_f)
end</
'''Testing:'''
<
[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321],
2)</
{{out}}
<pre>[1.0, 2.0, 3.0]</pre>
Line 1,815 ⟶ 2,338:
{{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)
Line 1,848 ⟶ 2,371:
polyRegression(xy)
}</
=={{header|Sidef}}==
{{trans|Ruby}}
<
var A = Matrix.build(x.len, degree+1, {|i,j|
x[i]**j
Line 1,871 ⟶ 2,394:
)
say coeff</
{{out}}
<pre>[1, 2, 3]</pre>
Line 1,877 ⟶ 2,400:
=={{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
Line 1,909 ⟶ 2,432:
|
_cons | 1 . . . . .
------------------------------------------------------------------------------</
=={{header|Swift}}==
{{trans|Kotlin}}
<syntaxhighlight lang="swift">
let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
Line 1,956 ⟶ 2,479:
polyRegression(x: x, y: y)
</syntaxhighlight>
{{out}}
Line 1,981 ⟶ 2,504:
<!-- This implementation from Emiliano Gavilan;
posted here with his explicit permission -->
<
proc build.matrix {xvec degree} {
Line 2,029 ⟶ 2,552:
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.
=={{header|TI-83 BASIC}}==
<syntaxhighlight lang="ti83b">DelVar X
seq(X,X,0,10) → L1
{1,6,17,34,57,86,121,162,209,262,321} → L2
QuadReg L1,L2</syntaxhighlight>
{{out}}
<pre>y=ax²+bx+c
a=3
b=2
c=1
</pre>
=={{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
Disp regeq(x)</
<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.
Line 2,057 ⟶ 2,594:
whereby the data can be passed as lists rather than arrays,
and all memory management is handled automatically.
<
#import nat
#import flo
(fit "n") ("x","y") = ..dgelsd\"y" (gang \/*pow float*x iota successor "n")* "x"</
test program:
<
y = <1.,6.,17.,34.,57.,86.,121.,162.,209.,262.,321.>
#cast %eL
example = fit2(x,y)</
{{out}}
<pre><3.000000e+00,2.000000e+00,1.000000e+00></pre>
Line 2,074 ⟶ 2,611:
=={{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
Line 2,103 ⟶ 2,640:
Debug.Print "Degrees of freedom:"; result(4, 2)
Debug.Print "Standard error of y estimate:"; result(3, 2)
End Sub</
<pre>coefficients : 1, 2, 3,
standard errors: 0, 0, 0,
Line 2,111 ⟶ 2,648:
Degrees of freedom: 8
Standard error of y estimate: 2,79482284961344E-14 </pre>
=={{header|Wren}}==
{{trans|REXX}}
{{libheader|Wren-math}}
{{libheader|Wren-seq}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./math" for Nums
import "./seq" for Lst
import "./fmt" for Fmt
var polynomialRegression = Fn.new { |x, y|
var xm = Nums.mean(x)
var ym = Nums.mean(y)
var x2m = Nums.mean(x.map { |e| e * e })
var x3m = Nums.mean(x.map { |e| e * e * e })
var x4m = Nums.mean(x.map { |e| e * e * e * e })
var z = Lst.zip(x, y)
var xym = Nums.mean(z.map { |p| p[0] * p[1] })
var x2ym = Nums.mean(z.map { |p| p[0] * p[0] * p[1] })
var sxx = x2m - xm * xm
var sxy = xym - xm * ym
var sxx2 = x3m - xm * x2m
var sx2x2 = x4m - x2m * x2m
var sx2y = x2ym - x2m * ym
var b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
var c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
var a = ym - b * xm - c * x2m
var abc = Fn.new { |xx| a + b * xx + c * xx * xx }
System.print("y = %(a) + %(b)x + %(c)x^2\n")
System.print(" Input Approximation")
System.print(" x y y1")
for (p in z) Fmt.print("$2d $3d $5.1f", p[0], p[1], abc.call(p[0]))
}
var x = List.filled(11, 0)
for (i in 1..10) x[i] = i
var y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]
polynomialRegression.call(x, y)</syntaxhighlight>
{{out}}
<pre>
y = 1 + 2x + 3x^2
Input Approximation
x y y1
0 1 1.0
1 6 6.0
2 17 17.0
3 34 34.0
4 57 57.0
5 86 86.0
6 121 121.0
7 162 162.0
8 209 209.0
9 262 262.0
10 321 321.0
</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);
Line 2,120 ⟶ 2,718:
v.format().println();
GSL.Helpers.polyString(v).println();
GSL.Helpers.polyEval(v,xs).format().println();</
{{out}}
<pre>
Line 2,132 ⟶ 2,730:
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)
.flatten().println();</
{{out}}<pre>L(1,2,3)</pre>
|