Polynomial regression: Difference between revisions
(add TI-89 BASIC) |
|||
(150 intermediate revisions by 69 users not shown) | |||
Line 1: | Line 1: | ||
{{task|Mathematical operations}} |
{{task|Mathematical operations|Matrices}} |
||
Find an approximating |
Find an approximating polynomial of known degree for a given data. |
||
Example: |
Example: |
||
Line 6: | Line 6: | ||
x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; |
x = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; |
||
y = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}; |
y = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}; |
||
The approximating |
The approximating polynomial is: |
||
3 x<sup>2</sup> + 2 x + 1 |
3 x<sup>2</sup> + 2 x + 1 |
||
Here, the |
Here, the polynomial's coefficients are (3, 2, 1). |
||
This task is intended as a subtask for [[Measure relative performance of sorting algorithms implementations]]. |
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}}== |
=={{header|Ada}}== |
||
<syntaxhighlight lang="ada">with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays; |
|||
<lang ada> |
|||
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays; |
|||
function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is |
function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is |
||
Line 26: | Line 80: | ||
end loop; |
end loop; |
||
return Solve (A * Transpose (A), A * Y); |
return Solve (A * Transpose (A), A * Y); |
||
end Fit; |
end Fit;</syntaxhighlight> |
||
</lang> |
|||
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. |
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=== |
===Example=== |
||
<syntaxhighlight lang="ada">with Fit; |
|||
<lang ada> |
|||
with Fit; |
|||
with Ada.Float_Text_IO; use Ada.Float_Text_IO; |
with Ada.Float_Text_IO; use Ada.Float_Text_IO; |
||
Line 45: | Line 97: | ||
Put (C (1), Aft => 3, Exp => 0); |
Put (C (1), Aft => 3, Exp => 0); |
||
Put (C (2), Aft => 3, Exp => 0); |
Put (C (2), Aft => 3, Exp => 0); |
||
end Fitting; |
end Fitting;</syntaxhighlight> |
||
{{out}} |
|||
</lang> |
|||
Sample output: |
|||
<pre> |
<pre> |
||
1.000 2.000 3.000 |
1.000 2.000 3.000 |
||
Line 55: | Line 106: | ||
{{trans|Ada}} |
{{trans|Ada}} |
||
{{works with|ALGOL 68|Standard - ''lu decomp'' and ''lu solve'' are from the [[:Category:Libgsl|GSL]] library}} |
{{works with|ALGOL 68|Standard - ''lu decomp'' and ''lu solve'' are from the [[:Category:Libgsl|GSL]] library}}<br> |
||
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}} <!-- with version's of a68g up to mk18 the GSL library has problems with VEC and MAT with LWB 0 , and transposes of non square MAT, hence we include source code of replacement OPerators here --> |
{{works with|ALGOL 68G|Any - tested with release mk15-0.8b.fc9.i386}} <!-- with version's of a68g up to mk18 the GSL library has problems with VEC and MAT with LWB 0 , and transposes of non square MAT, hence we include source code of replacement OPerators here --> |
||
<!-- {{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}} --> |
<!-- {{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}} --> |
||
< |
<syntaxhighlight lang="algol68">MODE FIELD = REAL; |
||
MODE |
MODE |
||
Line 171: | Line 223: | ||
); |
); |
||
print polynomial(d) |
print polynomial(d) |
||
END # fitting #</ |
END # fitting #</syntaxhighlight> |
||
{{out}} |
|||
Output: |
|||
<pre> |
<pre> |
||
3x**2+2x+1 |
3x**2+2x+1 |
||
1.0848x**2+10.3552x-0.6164 |
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> |
|||
=={{header|BBC BASIC}}== |
|||
{{works with|BBC BASIC for Windows}} |
|||
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! |
|||
<syntaxhighlight 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</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
</pre> |
||
=={{header|C}}== |
=={{header|C}}== |
||
{{libheader| |
{{libheader|GNU Scientific Library}} |
||
'''Include''' file (to make the code reusable easily) named <tt>polifitgsl.h</tt> |
'''Include''' file (to make the code reusable easily) named <tt>polifitgsl.h</tt> |
||
< |
<syntaxhighlight lang="c">#ifndef _POLIFITGSL_H |
||
#define _POLIFITGSL_H |
#define _POLIFITGSL_H |
||
#include <gsl/gsl_multifit.h> |
#include <gsl/gsl_multifit.h> |
||
Line 188: | Line 477: | ||
bool polynomialfit(int obs, int degree, |
bool polynomialfit(int obs, int degree, |
||
double *dx, double *dy, double *store); /* n, p */ |
double *dx, double *dy, double *store); /* n, p */ |
||
#endif</ |
#endif</syntaxhighlight> |
||
'''Implementation''' (the examples [http://www.gnu.org/software/gsl/manual/html_node/Fitting-Examples.html here] helped alot to code this quickly): |
'''Implementation''' (the examples [http://www.gnu.org/software/gsl/manual/html_node/Fitting-Examples.html here] helped alot to code this quickly): |
||
< |
<syntaxhighlight lang="c">#include "polifitgsl.h" |
||
bool polynomialfit(int obs, int degree, |
bool polynomialfit(int obs, int degree, |
||
Line 208: | Line 497: | ||
for(i=0; i < obs; i++) { |
for(i=0; i < obs; i++) { |
||
for(j=0; j < degree; j++) { |
|||
for(j=1; j < degree; j++) { |
|||
gsl_matrix_set(X, i, j, pow(dx[i], j)); |
gsl_matrix_set(X, i, j, pow(dx[i], j)); |
||
} |
} |
||
Line 231: | Line 519: | ||
return true; /* we do not "analyse" the result (cov matrix mainly) |
return true; /* we do not "analyse" the result (cov matrix mainly) |
||
to know if the fit is "good" */ |
to know if the fit is "good" */ |
||
}</ |
}</syntaxhighlight> |
||
'''Testing''': |
'''Testing''': |
||
< |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include "polifitgsl.h" |
#include "polifitgsl.h" |
||
Line 253: | Line 541: | ||
} |
} |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
|||
'''Output''' of the test: |
|||
<pre>1.000000 |
<pre>1.000000 |
||
2.000000 |
2.000000 |
||
3.000000</pre> |
3.000000</pre> |
||
=={{header|C sharp|C#}}== |
|||
{{libheader|Math.Net}} |
|||
<syntaxhighlight lang="csharp"> 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<double> 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(); |
|||
}</syntaxhighlight> |
|||
Example: |
|||
<syntaxhighlight 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], Polynomial.Evaluate(x[i], p), y[i] - Polynomial.Evaluate(x[i], p)); |
|||
Console.ReadKey(true); |
|||
}</syntaxhighlight> |
|||
=={{header|C++}}== |
|||
{{trans|Java}} |
|||
<syntaxhighlight lang="cpp">#include <algorithm> |
|||
#include <iostream> |
|||
#include <numeric> |
|||
#include <vector> |
|||
void polyRegression(const std::vector<int>& x, const std::vector<int>& y) { |
|||
int n = x.size(); |
|||
std::vector<int> r(n); |
|||
std::iota(r.begin(), r.end(), 0); |
|||
double xm = std::accumulate(x.begin(), x.end(), 0.0) / x.size(); |
|||
double ym = std::accumulate(y.begin(), y.end(), 0.0) / y.size(); |
|||
double x2m = std::transform_reduce(r.begin(), r.end(), 0.0, std::plus<double>{}, [](double a) {return a * a; }) / r.size(); |
|||
double x3m = std::transform_reduce(r.begin(), r.end(), 0.0, std::plus<double>{}, [](double a) {return a * a * a; }) / r.size(); |
|||
double x4m = std::transform_reduce(r.begin(), r.end(), 0.0, std::plus<double>{}, [](double a) {return a * a * a * a; }) / r.size(); |
|||
double xym = std::transform_reduce(x.begin(), x.end(), y.begin(), 0.0, std::plus<double>{}, std::multiplies<double>{}); |
|||
xym /= fmin(x.size(), y.size()); |
|||
double x2ym = std::transform_reduce(x.begin(), x.end(), y.begin(), 0.0, std::plus<double>{}, [](double a, double b) { return a * a * b; }); |
|||
x2ym /= fmin(x.size(), y.size()); |
|||
double sxx = x2m - xm * xm; |
|||
double sxy = xym - xm * ym; |
|||
double sxx2 = x3m - xm * x2m; |
|||
double sx2x2 = x4m - x2m * x2m; |
|||
double sx2y = x2ym - x2m * ym; |
|||
double b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2); |
|||
double c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2); |
|||
double a = ym - b * xm - c * x2m; |
|||
auto abc = [a, b, c](int xx) { |
|||
return a + b * xx + c * xx*xx; |
|||
}; |
|||
std::cout << "y = " << a << " + " << b << "x + " << c << "x^2" << std::endl; |
|||
std::cout << " Input Approximation" << std::endl; |
|||
std::cout << " x y y1" << std::endl; |
|||
auto xit = x.cbegin(); |
|||
auto xend = x.cend(); |
|||
auto yit = y.cbegin(); |
|||
auto yend = y.cend(); |
|||
while (xit != xend && yit != yend) { |
|||
printf("%2d %3d %5.1f\n", *xit, *yit, abc(*xit)); |
|||
xit = std::next(xit); |
|||
yit = std::next(yit); |
|||
} |
|||
} |
|||
int main() { |
|||
using namespace std; |
|||
vector<int> x(11); |
|||
iota(x.begin(), x.end(), 0); |
|||
vector<int> y{ 1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321 }; |
|||
polyRegression(x, y); |
|||
return 0; |
|||
}</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|Common Lisp}}== |
|||
Uses the routine (lsqr A b) from [[Multiple regression]] and (mtp A) from [[Matrix transposition]]. |
|||
<syntaxhighlight 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))))</syntaxhighlight> |
|||
Example: |
|||
<syntaxhighlight 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))</syntaxhighlight> |
|||
=={{header|D}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="d">import std.algorithm; |
|||
import std.range; |
|||
import std.stdio; |
|||
auto average(R)(R r) { |
|||
auto t = r.fold!("a+b", "a+1")(0, 0); |
|||
return cast(double) t[0] / t[1]; |
|||
} |
|||
void polyRegression(int[] x, int[] y) { |
|||
auto n = x.length; |
|||
auto r = iota(0, n).array; |
|||
auto xm = x.average(); |
|||
auto ym = y.average(); |
|||
auto x2m = r.map!"a*a".average(); |
|||
auto x3m = r.map!"a*a*a".average(); |
|||
auto x4m = r.map!"a*a*a*a".average(); |
|||
auto xym = x.zip(y).map!"a[0]*a[1]".average(); |
|||
auto x2ym = x.zip(y).map!"a[0]*a[0]*a[1]".average(); |
|||
auto sxx = x2m - xm * xm; |
|||
auto sxy = xym - xm * ym; |
|||
auto sxx2 = x3m - xm * x2m; |
|||
auto sx2x2 = x4m - x2m * x2m; |
|||
auto sx2y = x2ym - x2m * ym; |
|||
auto b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2); |
|||
auto c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2); |
|||
auto a = ym - b * xm - c * x2m; |
|||
real abc(int xx) { |
|||
return a + b * xx + c * xx * xx; |
|||
} |
|||
writeln("y = ", a, " + ", b, "x + ", c, "x^2"); |
|||
writeln(" Input Approximation"); |
|||
writeln(" x y y1"); |
|||
foreach (i; 0..n) { |
|||
writefln("%2d %3d %5.1f", x[i], y[i], abc(x[i])); |
|||
} |
|||
} |
|||
void main() { |
|||
auto x = iota(0, 11).array; |
|||
auto y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]; |
|||
polyRegression(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|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}}== |
=={{header|Fortran}}== |
||
{{libheader|LAPACK}} |
{{libheader|LAPACK}} |
||
< |
<syntaxhighlight lang="fortran">module fitting |
||
contains |
|||
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) |
|||
real(dp), dimension(:,:), allocatable :: X |
|||
real(dp), dimension(:,:), allocatable :: XT |
|||
real(dp), dimension(:,:), allocatable :: XTX |
|||
deallocate(ipiv) |
|||
integer :: i, j |
|||
deallocate(work) |
|||
deallocate(X) |
|||
integer :: n, lda, lwork |
|||
deallocate(XT) |
|||
integer :: info |
|||
deallocate(XTX) |
|||
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 function |
|||
</lang> |
|||
end module</syntaxhighlight> |
|||
===Example=== |
===Example=== |
||
<syntaxhighlight lang="fortran">program PolynomalFitting |
|||
<lang fortran> |
|||
use fitting |
|||
program PolynomalFitting |
|||
implicit none |
|||
use fitting |
|||
implicit none |
|||
! let us test it |
|||
integer, parameter :: degree = 2 |
|||
! let us test it |
|||
integer, parameter :: dp = selected_real_kind(15, 307) |
|||
integer :: i |
|||
real(dp), dimension(11) :: x = (/ (i,i=0,10) /) |
|||
integer :: i |
|||
real(dp), dimension(11) :: y = (/ 1, 6, 17, 34, & |
|||
57, 86, 121, 162, & |
|||
209, 262, 321 /) |
|||
real(dp), dimension(degree+1) :: a |
|||
209, 262, 321 /) |
|||
real(dp), dimension(degree+1) :: a |
|||
a = polyfit(x, y, degree) |
|||
a = polyfit(x, y, degree) |
|||
write (*, '(F9.4)') a |
|||
write (*, '(F9.4)') a |
|||
end program</syntaxhighlight> |
|||
</lang> |
|||
Output (lower powers first, so this seems the opposite of the Python output): |
|||
{{out}} (lower powers first, so this seems the opposite of the Python output): |
|||
<pre> |
<pre> |
||
1.0000 |
1.0000 |
||
Line 358: | Line 894: | ||
3.0000 |
3.0000 |
||
</pre> |
</pre> |
||
=={{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 matrix |
|||
Dim As Double element(Any,Any) |
|||
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 |
|||
Next z |
|||
Function = ans |
|||
End Function |
|||
Function matrix.inverse() As matrix |
|||
Var ub1=Ubound(this.element,1),ub2=Ubound(this.element,2) |
|||
Dim As matrix ans |
|||
Dim As vector temp,null_ |
|||
Redim temp.element(1 To ub1):Redim null_.element(1 To ub1) |
|||
Redim ans.element(1 To ub1,1 To ub2) |
|||
For a As Integer=1 To ub1 |
|||
temp=null_ |
|||
temp.element(a)=1 |
|||
temp=GaussJordan(temp) |
|||
For b As Integer=1 To ub1 |
|||
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}}== |
|||
<syntaxhighlight 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 ]</syntaxhighlight> |
|||
=={{header|gnuplot}}== |
=={{header|gnuplot}}== |
||
< |
<syntaxhighlight lang="gnuplot"># The polynomial approximation |
||
f(x) = a*x**2 + b*x + c |
f(x) = a*x**2 + b*x + c |
||
Line 384: | Line 1,127: | ||
e |
e |
||
print sprintf("\n --- \n Polynomial fit: %.4f x^2 + %.4f x + %.4f\n", a, b, c)</ |
print sprintf("\n --- \n Polynomial fit: %.4f x^2 + %.4f x + %.4f\n", a, b, c)</syntaxhighlight> |
||
=={{header|Go}}== |
|||
===Library gonum/matrix=== |
|||
<syntaxhighlight lang="go">package main |
|||
import ( |
|||
"fmt" |
|||
"log" |
|||
"gonum.org/v1/gonum/mat" |
|||
) |
|||
func main() { |
|||
var ( |
|||
x = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10} |
|||
y = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321} |
|||
degree = 2 |
|||
a = Vandermonde(x, degree+1) |
|||
b = mat.NewDense(len(y), 1, y) |
|||
c = mat.NewDense(degree+1, 1, nil) |
|||
) |
|||
var qr mat.QR |
|||
qr.Factorize(a) |
|||
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, d int) *mat.Dense { |
|||
x := mat.NewDense(len(a), d, nil) |
|||
for i := range a { |
|||
for j, p := 0, 1.0; j < d; j, p = j+1, p*a[i] { |
|||
x.Set(i, j, p) |
|||
} |
|||
} |
|||
return x |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
⎡1.000⎤ |
|||
⎢2.000⎥ |
|||
⎣3.000⎦ |
|||
</pre> |
|||
===Library go.matrix=== |
|||
Least squares solution using QR decomposition and package [http://github.com/skelterjohn/go.matrix go.matrix]. |
|||
<syntaxhighlight lang="go">package main |
|||
import ( |
|||
"fmt" |
|||
"github.com/skelterjohn/go.matrix" |
|||
) |
|||
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) |
|||
}</syntaxhighlight> |
|||
{{out}} (lowest order coefficient first) |
|||
<pre> |
|||
[0.9999999999999758 2.000000000000015 2.999999999999999] |
|||
</pre> |
|||
=={{header|Haskell}}== |
|||
Uses module Matrix.LU from [http://hackage.haskell.org/package/dsp hackageDB DSP] |
|||
<syntaxhighlight 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</syntaxhighlight> |
|||
{{out}} in GHCi: |
|||
<syntaxhighlight lang="haskell">*Main> polyfit 3 [1,6,17,34,57,86,121,162,209,262,321] |
|||
[1.0,2.0,3.0]</syntaxhighlight> |
|||
=={{header|HicEst}}== |
|||
<syntaxhighlight 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</syntaxhighlight> |
|||
{{out}} |
|||
<pre>SOLVE performs a (nonlinear) least-square fit (Levenberg-Marquardt): |
|||
p(1)=2.997135145; p(2)=2.011348347; p(3)=0.9906627242; iterations=19;</pre> |
|||
=={{header|Hy}}== |
|||
<syntaxhighlight lang="lisp">(import [numpy [polyfit]]) |
|||
(setv x (range 11)) |
|||
(setv y [1 6 17 34 57 86 121 162 209 262 321]) |
|||
(print (polyfit x y 2))</syntaxhighlight> |
|||
=={{header|J}}== |
|||
<syntaxhighlight lang="j"> Y=:1 6 17 34 57 86 121 162 209 262 321 |
|||
(%. ^/~@x:@i.@#) Y |
|||
1 2 3 0 0 0 0 0 0 0 0</syntaxhighlight> |
|||
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 (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: |
|||
<syntaxhighlight lang="j"> Y %. (i.3) ^/~ i.#Y |
|||
1 2 3</syntaxhighlight> |
|||
(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|Java}}== |
|||
{{trans|D}} |
|||
{{works with|Java|8}} |
|||
<syntaxhighlight lang="java">import java.util.Arrays; |
|||
import java.util.function.IntToDoubleFunction; |
|||
import java.util.stream.IntStream; |
|||
public class PolynomialRegression { |
|||
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(x).map(a -> a * a).average().orElse(Double.NaN); |
|||
double x3m = Arrays.stream(x).map(a -> a * a * a).average().orElse(Double.NaN); |
|||
double x4m = Arrays.stream(x).map(a -> a * a * a * a).average().orElse(Double.NaN); |
|||
double xym = 0.0; |
|||
for (int i = 0; i < x.length && i < y.length; ++i) { |
|||
xym += x[i] * y[i]; |
|||
} |
|||
xym /= Math.min(x.length, y.length); |
|||
double x2ym = 0.0; |
|||
for (int i = 0; i < x.length && i < y.length; ++i) { |
|||
x2ym += x[i] * x[i] * y[i]; |
|||
} |
|||
x2ym /= Math.min(x.length, y.length); |
|||
double sxx = x2m - xm * xm; |
|||
double sxy = xym - xm * ym; |
|||
double sxx2 = x3m - xm * x2m; |
|||
double sx2x2 = x4m - x2m * x2m; |
|||
double sx2y = x2ym - x2m * ym; |
|||
double b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2); |
|||
double c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2); |
|||
double a = ym - b * xm - c * x2m; |
|||
IntToDoubleFunction abc = (int xx) -> a + b * xx + c * xx * xx; |
|||
System.out.println("y = " + a + " + " + b + "x + " + c + "x^2"); |
|||
System.out.println(" Input Approximation"); |
|||
System.out.println(" x y y1"); |
|||
for (int i = 0; i < n; ++i) { |
|||
System.out.printf("%2d %3d %5.1f\n", x[i], y[i], abc.applyAsDouble(x[i])); |
|||
} |
|||
} |
|||
public static void main(String[] args) { |
|||
int[] x = IntStream.range(0, 11).toArray(); |
|||
int[] y = new int[]{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}; |
|||
polyRegression(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|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}}== |
|||
{{works with|Julia|0.6}} |
|||
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): |
|||
<syntaxhighlight lang="julia">polyfit(x::Vector, y::Vector, deg::Int) = collect(v ^ p for v in x, p in 0:deg) \ y |
|||
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)</syntaxhighlight> |
|||
{{out}} |
|||
<pre>polyfit(x, y, 2) = [1.0, 2.0, 3.0]</pre> |
|||
=={{header|Kotlin}}== |
|||
{{trans|REXX}} |
|||
<syntaxhighlight lang="scala">// version 1.1.51 |
|||
fun polyRegression(x: IntArray, y: IntArray) { |
|||
val xm = x.average() |
|||
val ym = y.average() |
|||
val x2m = x.map { it * it }.average() |
|||
val x3m = x.map { it * it * it }.average() |
|||
val x4m = x.map { it * it * it * it }.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() |
|||
val sxx = x2m - xm * xm |
|||
val sxy = xym - xm * ym |
|||
val sxx2 = x3m - xm * x2m |
|||
val sx2x2 = x4m - x2m * x2m |
|||
val sx2y = x2ym - x2m * ym |
|||
val b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) |
|||
val c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) |
|||
val a = ym - b * xm - c * x2m |
|||
fun abc(xx: Int) = a + b * xx + c * xx * xx |
|||
println("y = $a + ${b}x + ${c}x^2\n") |
|||
println(" Input Approximation") |
|||
println(" x y y1") |
|||
for ((xi, yi) in x zip y) { |
|||
System.out.printf("%2d %3d %5.1f\n", xi, yi, abc(xi)) |
|||
} |
|||
} |
|||
fun main() { |
|||
val x = IntArray(11) { it } |
|||
val y = intArrayOf(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321) |
|||
polyRegression(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|Lua}}== |
|||
{{trans|Modula-2}} |
|||
<syntaxhighlight lang="lua">function eval(a,b,c,x) |
|||
return a + (b + c * x) * x |
|||
end |
|||
function regression(xa,ya) |
|||
local n = #xa |
|||
local xm = 0.0 |
|||
local ym = 0.0 |
|||
local x2m = 0.0 |
|||
local x3m = 0.0 |
|||
local x4m = 0.0 |
|||
local xym = 0.0 |
|||
local x2ym = 0.0 |
|||
for i=1,n do |
|||
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] |
|||
end |
|||
xm = xm / n |
|||
ym = ym / n |
|||
x2m = x2m / n |
|||
x3m = x3m / n |
|||
x4m = x4m / n |
|||
xym = xym / n |
|||
x2ym = x2ym / n |
|||
local sxx = x2m - xm * xm |
|||
local sxy = xym - xm * ym |
|||
local sxx2 = x3m - xm * x2m |
|||
local sx2x2 = x4m - x2m * x2m |
|||
local sx2y = x2ym - x2m * ym |
|||
local b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) |
|||
local c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) |
|||
local a = ym - b * xm - c * x2m |
|||
print("y = "..a.." + "..b.."x + "..c.."x^2") |
|||
for i=1,n do |
|||
print(string.format("%2d %3d %3d", xa[i], ya[i], eval(a, b, c, xa[i]))) |
|||
end |
|||
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} |
|||
regression(xa, ya)</syntaxhighlight> |
|||
{{out}} |
|||
<pre>y = 1 + 2x + 3x^2 |
|||
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|Maple}}== |
|||
<syntaxhighlight lang="maple">with(CurveFitting); |
|||
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. |
|||
<syntaxhighlight lang="mathematica">Simplify@InterpolatingPolynomial[{{0, 1}, {1, 6}, {2, 17}, {3, 34}, {4, 57}, {5, 86}, {6, 121}, {7, 162}, {8, 209}, {9, 262}, {10, 321}}, x]</syntaxhighlight> |
|||
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> |
|||
=={{header|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 <math>point_{i} = (x_{i},y_{i})</math> 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. |
|||
<syntaxhighlight 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</syntaxhighlight> |
|||
=={{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</syntaxhighlight> |
|||
''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 * - |
|||
ПA ИП4 ИП3 * ИП2 ИП5 * - П3 ИП4 |
|||
ИП8 * ИП2 ИП6 * - П8 ИП4 ИПB * |
|||
ИП2 ИП9 * - ИПD * ИП3 ИПA * - |
|||
ИПD ИП8 * ИП7 ИП3 * - / ПB ИПA |
|||
ИПB ИП7 * - ИПD / ПA ИП9 ИПB ИП6 |
|||
* - ИПA ИП5 * - ИП4 / П9 С/П</syntaxhighlight> |
|||
''Result'': Р9 = a<sub>0</sub>, РA = a<sub>1</sub>, РB = a<sub>2</sub>. |
|||
=={{header|Modula-2}}== |
|||
<syntaxhighlight lang="modula2">MODULE PolynomialRegression; |
|||
FROM FormatString IMPORT FormatString; |
|||
FROM RealStr IMPORT RealToStr; |
|||
FROM Terminal IMPORT WriteString,WriteLn,ReadChar; |
|||
PROCEDURE Eval(a,b,c,x : REAL) : REAL; |
|||
BEGIN |
|||
RETURN a + b*x + c*x*x; |
|||
END Eval; |
|||
PROCEDURE Regression(x,y : ARRAY OF INTEGER); |
|||
VAR |
|||
n,i : INTEGER; |
|||
xm,x2m,x3m,x4m : REAL; |
|||
ym : REAL; |
|||
xym,x2ym : REAL; |
|||
sxx,sxy,sxx2,sx2x2,sx2y : REAL; |
|||
a,b,c : REAL; |
|||
buf : ARRAY[0..63] OF CHAR; |
|||
BEGIN |
|||
n := SIZE(x)/SIZE(INTEGER); |
|||
xm := 0.0; |
|||
ym := 0.0; |
|||
x2m := 0.0; |
|||
x3m := 0.0; |
|||
x4m := 0.0; |
|||
xym := 0.0; |
|||
x2ym := 0.0; |
|||
FOR i:=0 TO n-1 DO |
|||
xm := xm + FLOAT(x[i]); |
|||
ym := ym + FLOAT(y[i]); |
|||
x2m := x2m + FLOAT(x[i]) * FLOAT(x[i]); |
|||
x3m := x3m + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]); |
|||
x4m := x4m + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]); |
|||
xym := xym + FLOAT(x[i]) * FLOAT(y[i]); |
|||
x2ym := x2ym + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(y[i]); |
|||
END; |
|||
xm := xm / FLOAT(n); |
|||
ym := ym / FLOAT(n); |
|||
x2m := x2m / FLOAT(n); |
|||
x3m := x3m / FLOAT(n); |
|||
x4m := x4m / FLOAT(n); |
|||
xym := xym / FLOAT(n); |
|||
x2ym := x2ym / FLOAT(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; |
|||
WriteString("y = "); |
|||
RealToStr(a, buf); |
|||
WriteString(buf); |
|||
WriteString(" + "); |
|||
RealToStr(b, buf); |
|||
WriteString(buf); |
|||
WriteString("x + "); |
|||
RealToStr(c, buf); |
|||
WriteString(buf); |
|||
WriteString("x^2"); |
|||
WriteLn; |
|||
FOR i:=0 TO n-1 DO |
|||
FormatString("%2i %3i ", buf, x[i], y[i]); |
|||
WriteString(buf); |
|||
RealToStr(Eval(a,b,c,FLOAT(x[i])), buf); |
|||
WriteString(buf); |
|||
WriteLn; |
|||
END; |
|||
END Regression; |
|||
TYPE R = ARRAY[0..10] OF INTEGER; |
|||
VAR |
|||
x,y : R; |
|||
BEGIN |
|||
x := R{0,1,2,3,4,5,6,7,8,9,10}; |
|||
y := R{1,6,17,34,57,86,121,162,209,262,321}; |
|||
Regression(x,y); |
|||
ReadChar; |
|||
END PolynomialRegression.</syntaxhighlight> |
|||
=={{header|Nim}}== |
|||
{{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}}== |
=={{header|Octave}}== |
||
< |
<syntaxhighlight lang="octave">x = [0:10]; |
||
y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]; |
y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]; |
||
coeffs = polyfit(x, y, 2)</ |
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.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 |
|||
x y f(x) |
|||
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|PowerShell}}== |
|||
<syntaxhighlight lang="powershell"> |
|||
function qr([double[][]]$A) { |
|||
$m,$n = $A.count, $A[0].count |
|||
$pm,$pn = ($m-1), ($n-1) |
|||
[double[][]]$Q = 0..($m-1) | foreach{$row = @(0) * $m; $row[$_] = 1; ,$row} |
|||
[double[][]]$R = $A | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})} |
|||
foreach ($h in 0..$pn) { |
|||
[double[]]$u = $R[$h..$pm] | foreach{$_[$h]} |
|||
[double]$nu = $u | foreach {[double]$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)} |
|||
$u[0] -= if ($u[0] -lt 1) {$nu} else {-$nu} |
|||
[double]$nu = $u | foreach {$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)} |
|||
[double[]]$u = $u | foreach { $_/$nu} |
|||
[double[][]]$v = 0..($u.Count - 1) | foreach{$i = $_; ,($u | foreach{2*$u[$i]*$_})} |
|||
[double[][]]$CR = $R | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})} |
|||
[double[][]]$CQ = $Q | foreach{$row = $_; ,@(0..$pm | foreach{$row[$_]})} |
|||
foreach ($i in $h..$pm) { |
|||
foreach ($j in $h..$pn) { |
|||
$R[$i][$j] -= $h..$pm | foreach {[double]$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CR[$_][$j]} {$sum} |
|||
} |
|||
} |
|||
if (0 -eq $h) { |
|||
foreach ($i in $h..$pm) { |
|||
foreach ($j in $h..$pm) { |
|||
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i][$_]*$CQ[$_][$j]} {$sum} |
|||
} |
|||
} |
|||
} else { |
|||
$p = $h-1 |
|||
foreach ($i in $h..$pm) { |
|||
foreach ($j in 0..$p) { |
|||
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum} |
|||
} |
|||
foreach ($j in $h..$pm) { |
|||
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum} |
|||
} |
|||
} |
|||
} |
|||
} |
|||
foreach ($i in 0..$pm) { |
|||
foreach ($j in $i..$pm) {$Q[$i][$j],$Q[$j][$i] = $Q[$j][$i],$Q[$i][$j]} |
|||
} |
|||
[PSCustomObject]@{"Q" = $Q; "R" = $R} |
|||
} |
|||
function leastsquares([Double[][]]$A,[Double[]]$y) { |
|||
$QR = qr $A |
|||
[Double[][]]$Q = $QR.Q |
|||
[Double[][]]$R = $QR.R |
|||
$m,$n = $A.count, $A[0].count |
|||
[Double[]]$z = foreach ($j in 0..($m-1)) { |
|||
0..($m-1) | foreach {$sum = 0} {$sum += $Q[$_][$j]*$y[$_]} {$sum} |
|||
} |
|||
[Double[]]$x = @(0)*$n |
|||
for ($i = $n-1; $i -ge 0; $i--) { |
|||
for ($j = $i+1; $j -lt $n; $j++) { |
|||
$z[$i] -= $x[$j]*$R[$i][$j] |
|||
} |
|||
$x[$i] = $z[$i]/$R[$i][$i] |
|||
} |
|||
$x |
|||
} |
|||
function polyfit([Double[]]$x,[Double[]]$y,$n) { |
|||
$m = $x.Count |
|||
[Double[][]]$A = 0..($m-1) | foreach{$row = @(1) * ($n+1); ,$row} |
|||
for ($i = 0; $i -lt $m; $i++) { |
|||
for ($j = $n-1; 0 -le $j; $j--) { |
|||
$A[$i][$j] = $A[$i][$j+1]*$x[$i] |
|||
} |
|||
} |
|||
leastsquares $A $y |
|||
} |
|||
function show($m) {$m | foreach {write-host "$_"}} |
|||
$A = @(@(12,-51,4), @(6,167,-68), @(-4,24,-41)) |
|||
$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^2 X constant" |
|||
"$(polyfit $x $y 2)" |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
polyfit |
|||
X^2 X constant |
|||
3 1.99999999999998 1.00000000000005 |
|||
</pre> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
||
{{libheader| |
{{libheader|NumPy}} |
||
<syntaxhighlight lang="python">>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] |
|||
<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] |
>>> y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321] |
||
>>> coeffs = numpy.polyfit(x,y,deg=2) |
>>> coeffs = numpy.polyfit(x,y,deg=2) |
||
>>> coeffs |
>>> coeffs |
||
array([ 3., 2., 1.]) |
array([ 3., 2., 1.])</syntaxhighlight> |
||
</lang> |
|||
Substitute back received coefficients. |
Substitute back received coefficients. |
||
<syntaxhighlight lang="python">>>> yf = numpy.polyval(numpy.poly1d(coeffs), x) |
|||
<lang python> |
|||
>>> yf = numpy.polyval(numpy.poly1d(coeffs), x) |
|||
>>> yf |
>>> yf |
||
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.]) |
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.])</syntaxhighlight> |
||
Find max absolute error: |
|||
</lang> |
|||
<syntaxhighlight lang="python">>>> '%.1g' % max(y-yf) |
|||
Find max absolute error. |
|||
'1e-013'</syntaxhighlight> |
|||
<lang python> |
|||
>>> '%.1g' % max(y-yf) |
|||
'1e-013' |
|||
</lang> |
|||
===Example=== |
===Example=== |
||
For input arrays `x' and `y': |
For input arrays `x' and `y': |
||
<syntaxhighlight lang="python">>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] |
|||
<lang python> |
|||
>>> |
>>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0]</syntaxhighlight> |
||
>>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0] |
|||
</lang> |
|||
<syntaxhighlight lang="python">>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N') |
|||
<lang python> |
|||
>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N') |
|||
>>> print p |
>>> print p |
||
2 |
2 |
||
1.085 N + 10.36 N - 0.6164 |
1.085 N + 10.36 N - 0.6164</syntaxhighlight> |
||
Thus we confirm once more that for already sorted sequences |
|||
</lang> |
|||
the considered quick sort implementation has |
|||
Thus we confirm once more that for already sorted sequences the considered quick sort implementation has quadratic dependence on sequence length (see [[Query Performance|'''Example''' section for Python language on ''Query Performance'' page]]). |
|||
quadratic dependence on sequence length |
|||
(see [[Query Performance|'''Example''' section for Python language |
|||
on ''Query Performance'' page]]). |
|||
=={{header|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: |
|||
<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)))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
(Intercept) x I(x^2) |
|||
1 2 3 |
|||
</pre> |
|||
Alternately, use poly: |
|||
<syntaxhighlight lang="r">coef(lm(y ~ poly(x, 2, raw=T)))</syntaxhighlight>{{out}} |
|||
<pre> (Intercept) poly(x, 2, raw = T)1 poly(x, 2, raw = T)2 |
|||
1 2 3</pre> |
|||
=={{header|Racket}}== |
|||
<syntaxhighlight 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))))) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
[[File:polyreg-racket.png]] |
|||
=={{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>; |
|||
constant @y = <1 6 17 34 57 86 121 162 209 262 321>; |
|||
constant $x0 = [+] @e[^@x1]; |
|||
constant $x1 = [+] @x1 Z* @e; |
|||
constant $x2 = [+] @x1 »**» 2 Z* @e; |
|||
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 |
|||
2 |
|||
3 |
|||
</pre> |
|||
=={{header|REXX}}== |
|||
<syntaxhighlight lang="rexx">/* REXX --------------------------------------------------------------- |
|||
* Implementation of http://keisan.casio.com/exec/system/14059932254941 |
|||
*--------------------------------------------------------------------*/ |
|||
xl='0 1 2 3 4 5 6 7 8 9 10' |
|||
yl='1 6 17 34 57 86 121 162 209 262 321' |
|||
n=11 |
|||
Do i=1 To n |
|||
Parse Var xl x.i xl |
|||
Parse Var yl y.i yl |
|||
End |
|||
xm=0 |
|||
ym=0 |
|||
x2m=0 |
|||
x3m=0 |
|||
x4m=0 |
|||
xym=0 |
|||
x2ym=0 |
|||
Do i=1 To n |
|||
xm=xm+x.i |
|||
ym=ym+y.i |
|||
x2m=x2m+x.i**2 |
|||
x3m=x3m+x.i**3 |
|||
x4m=x4m+x.i**4 |
|||
xym=xym+x.i*y.i |
|||
x2ym=x2ym+(x.i**2)*y.i |
|||
End |
|||
xm =xm /n |
|||
ym =ym /n |
|||
x2m=x2m/n |
|||
x3m=x3m/n |
|||
x4m=x4m/n |
|||
xym=xym/n |
|||
x2ym=x2ym/n |
|||
Sxx=x2m-xm**2 |
|||
Sxy=xym-xm*ym |
|||
Sxx2=x3m-xm*x2m |
|||
Sx2x2=x4m-x2m**2 |
|||
Sx2y=x2ym-x2m*ym |
|||
B=(Sxy*Sx2x2-Sx2y*Sxx2)/(Sxx*Sx2x2-Sxx2**2) |
|||
C=(Sx2y*Sxx-Sxy*Sxx2)/(Sxx*Sx2x2-Sxx2**2) |
|||
A=ym-B*xm-C*x2m |
|||
Say 'y='a'+'||b'*x+'c'*x**2' |
|||
Say ' Input "Approximation"' |
|||
Say ' x y y1' |
|||
Do i=1 To 11 |
|||
Say right(x.i,2) right(y.i,3) format(fun(x.i),5,3) |
|||
End |
|||
Exit |
|||
fun: |
|||
Parse Arg x |
|||
Return a+b*x+c*x**2 </syntaxhighlight> |
|||
{{out}} |
|||
<pre>y=1+2*x+3*x**2 |
|||
Input "Approximation" |
|||
x y y1 |
|||
0 1 1.000 |
|||
1 6 6.000 |
|||
2 17 17.000 |
|||
3 34 34.000 |
|||
4 57 57.000 |
|||
5 86 86.000 |
|||
6 121 121.000 |
|||
7 162 162.000 |
|||
8 209 209.000 |
|||
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}}== |
|||
<syntaxhighlight lang="ruby">require 'matrix' |
|||
def regress x, y, degree |
|||
x_data = x.map { |xi| (0..degree).map { |pow| (xi**pow).to_r } } |
|||
mx = Matrix[*x_data] |
|||
my = Matrix.column_vector(y) |
|||
((mx.t * mx).inv * mx.t * my).transpose.to_a[0].map(&:to_f) |
|||
end</syntaxhighlight> |
|||
'''Testing:''' |
|||
<syntaxhighlight lang="ruby">p regress([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], |
|||
[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321], |
|||
2)</syntaxhighlight> |
|||
{{out}} |
|||
<pre>[1.0, 2.0, 3.0]</pre> |
|||
=={{header|Scala}}== |
|||
{{Out}}See it yourself by running in your browser [https://scastie.scala-lang.org/NklZH2LlScCpfsN4NSfFvA Scastie (remote JVM)]. |
|||
{{libheader|Scala Math Polynomial}} |
|||
{{libheader|Scastie qualified}} |
|||
{{works with|Scala|2.13}} |
|||
<syntaxhighlight lang="scala">object PolynomialRegression extends App { |
|||
private def xy = Seq(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321).zipWithIndex.map(_.swap) |
|||
private def polyRegression(xy: Seq[(Int, Int)]): Unit = { |
|||
val r = xy.indices |
|||
def average[U](ts: Iterable[U])(implicit num: Numeric[U]) = num.toDouble(ts.sum) / ts.size |
|||
def x3m: Double = average(r.map(a => a * a * a)) |
|||
def x4m: Double = average(r.map(a => a * a * a * a)) |
|||
def x2ym = xy.reduce((a, x) => (a._1 + x._1 * x._1 * x._2, 0))._1.toDouble / xy.size |
|||
def xym = xy.reduce((a, x) => (a._1 + x._1 * x._2, 0))._1.toDouble / xy.size |
|||
val x2m: Double = average(r.map(a => a * a)) |
|||
val (xm, ym) = (average(xy.map(_._1)), average(xy.map(_._2))) |
|||
val (sxx, sxy) = (x2m - xm * xm, xym - xm * ym) |
|||
val sxx2: Double = x3m - xm * x2m |
|||
val sx2x2: Double = x4m - x2m * x2m |
|||
val sx2y: Double = x2ym - x2m * ym |
|||
val c: Double = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) |
|||
val b: Double = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2) |
|||
val a: Double = ym - b * xm - c * x2m |
|||
def abc(xx: Int) = a + b * xx + c * xx * xx |
|||
println(s"y = $a + ${b}x + ${c}x^2") |
|||
println(" Input Approximation") |
|||
println(" x y y1") |
|||
xy.foreach {el => println(f"${el._1}%2d ${el._2}%3d ${abc(el._1)}%5.1f")} |
|||
} |
|||
polyRegression(xy) |
|||
}</syntaxhighlight> |
|||
=={{header|Sidef}}== |
|||
{{trans|Ruby}} |
|||
<syntaxhighlight lang="ruby">func regress(x, y, degree) { |
|||
var A = Matrix.build(x.len, degree+1, {|i,j| |
|||
x[i]**j |
|||
}) |
|||
var B = Matrix.column_vector(y...) |
|||
((A.transpose * A)**(-1) * A.transpose * B).transpose[0] |
|||
} |
|||
func poly(x) { |
|||
3*x**2 + 2*x + 1 |
|||
} |
|||
var coeff = regress( |
|||
10.of { _ }, |
|||
10.of { poly(_) }, |
|||
2 |
|||
) |
|||
say coeff</syntaxhighlight> |
|||
{{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. |
|||
<syntaxhighlight lang="stata">. clear |
|||
. input x y |
|||
0 1 |
|||
1 6 |
|||
2 17 |
|||
3 34 |
|||
4 57 |
|||
5 86 |
|||
6 121 |
|||
7 162 |
|||
8 209 |
|||
9 262 |
|||
10 321 |
|||
end |
|||
. regress y c.x##c.x |
|||
Source | SS df MS Number of obs = 11 |
|||
-------------+---------------------------------- F(2, 8) = . |
|||
Model | 120362 2 60181 Prob > F = . |
|||
Residual | 0 8 0 R-squared = 1.0000 |
|||
-------------+---------------------------------- Adj R-squared = 1.0000 |
|||
Total | 120362 10 12036.2 Root MSE = 0 |
|||
------------------------------------------------------------------------------ |
|||
y | Coef. Std. Err. t P>|t| [95% Conf. Interval] |
|||
-------------+---------------------------------------------------------------- |
|||
x | 2 . . . . . |
|||
| |
|||
c.x#c.x | 3 . . . . . |
|||
| |
|||
_cons | 1 . . . . . |
|||
------------------------------------------------------------------------------</syntaxhighlight> |
|||
=={{header|Swift}}== |
|||
{{trans|Kotlin}} |
|||
<syntaxhighlight lang="swift"> |
|||
let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] |
|||
let y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321] |
|||
func average(_ input: [Int]) -> Int { |
|||
return input.reduce(0, +) / input.count |
|||
} |
|||
func polyRegression(x: [Int], y: [Int]) { |
|||
let xm = average(x) |
|||
let ym = average(y) |
|||
let x2m = average(x.map { $0 * $0 }) |
|||
let x3m = average(x.map { $0 * $0 * $0 }) |
|||
let x4m = average(x.map { $0 * $0 * $0 * $0 }) |
|||
let xym = average(zip(x,y).map { $0 * $1 }) |
|||
let x2ym = average(zip(x,y).map { $0 * $0 * $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(xx: Int) -> Int { |
|||
return (a + b * xx) + (c * xx * xx) |
|||
} |
|||
print("y = \(a) + \(b)x + \(c)x^2\n") |
|||
print(" Input Approximation") |
|||
print(" x y y1") |
|||
for i in 0 ..< x.count { |
|||
let result = Double(abc(xx: i)) |
|||
print(String(format: "%2d %3d %5.1f", x[i], y[i], result)) |
|||
} |
|||
} |
|||
polyRegression(x: x, y: 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|Tcl}}== |
=={{header|Tcl}}== |
||
{{tcllib|math::linearalgebra}} |
|||
{{libheader|tcllib}}(which includes a package for performing linear algebra operations) |
|||
<!-- This implementation from Emiliano Gavilan; |
<!-- This implementation from Emiliano Gavilan; |
||
posted here with his explicit permission --> |
|||
<lang tcl>package require math::linearalgebra |
|||
<syntaxhighlight lang="tcl">package require math::linearalgebra |
|||
proc build.matrix {xvec degree} { |
proc build.matrix {xvec degree} { |
||
Line 480: | Line 2,552: | ||
set coeffs [math::linearalgebra::solveGauss $A $b] |
set coeffs [math::linearalgebra::solveGauss $A $b] |
||
# show results |
# show results |
||
puts $coeffs</ |
puts $coeffs</syntaxhighlight> |
||
This will print: |
This will print: |
||
1.0000000000000207 1.9999999999999958 3.0 |
1.0000000000000207 1.9999999999999958 3.0 |
||
which is a close approximation to the correct solution. |
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}}== |
=={{header|TI-89 BASIC}}== |
||
<syntaxhighlight lang="ti89b">DelVar x |
|||
<lang ti-89>seq(x,x,0,10)→xs |
|||
seq(x,x,0,10) → xs |
|||
{1,6,17,34,57,86,121,162,209,262,321}→ys |
|||
{1,6,17,34,57,86,121,162,209,262,321} → ys |
|||
QuadReg xs,ys |
QuadReg xs,ys |
||
Disp regeq(x)</ |
Disp regeq(x)</syntaxhighlight> |
||
<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. |
|||
Output: <code>3.·x<sup>2</sup> + 2.·x + 1.</code> |
|||
<code>QuadReg</code>, "quadratic regression", does the fit and stores the details in a number of standard variables, including <var>regeq</var>, which receives the fitted quadratic (polynomial) function itself. |
|||
We then apply that function to the (undefined as ensured by <code>DelVar</code>) variable x to obtain the expression in terms of x, and display it. |
|||
{{out}} |
|||
<code>3.·x<sup>2</sup> + 2.·x + 1.</code> |
|||
=={{header|Ursala}}== |
|||
{{libheader|LAPACK}} |
|||
The fit function defined below returns the coefficients |
|||
of an nth-degree polynomial in order of descending degree |
|||
fitting the lists of inputs x and outputs y. |
|||
The real work is done by the dgelsd function from the lapack library. |
|||
Ursala provides a simplified interface to this library |
|||
whereby the data can be passed as lists rather than arrays, |
|||
and all memory management is handled automatically. |
|||
<syntaxhighlight lang="ursala">#import std |
|||
#import nat |
|||
#import flo |
|||
(fit "n") ("x","y") = ..dgelsd\"y" (gang \/*pow float*x iota successor "n")* "x"</syntaxhighlight> |
|||
test program: |
|||
<syntaxhighlight 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)</syntaxhighlight> |
|||
{{out}} |
|||
<pre><3.000000e+00,2.000000e+00,1.000000e+00></pre> |
|||
=={{header|VBA}}== |
|||
Excel VBA has built in capability for line estimation. |
|||
<syntaxhighlight lang="vb">Option Base 1 |
|||
Private Function polynomial_regression(y As Variant, x As Variant, degree As Integer) As Variant |
|||
Dim a() As Double |
|||
ReDim a(UBound(x), 2) |
|||
For i = 1 To UBound(x) |
|||
For j = 1 To degree |
|||
a(i, j) = x(i) ^ j |
|||
Next j |
|||
Next i |
|||
polynomial_regression = WorksheetFunction.LinEst(WorksheetFunction.Transpose(y), a, True, True) |
|||
End Function |
|||
Public Sub main() |
|||
x = [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}] |
|||
y = [{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}] |
|||
result = polynomial_regression(y, x, 2) |
|||
Debug.Print "coefficients : "; |
|||
For i = UBound(result, 2) To 1 Step -1 |
|||
Debug.Print Format(result(1, i), "0.#####"), |
|||
Next i |
|||
Debug.Print |
|||
Debug.Print "standard errors: "; |
|||
For i = UBound(result, 2) To 1 Step -1 |
|||
Debug.Print Format(result(2, i), "0.#####"), |
|||
Next i |
|||
Debug.Print vbCrLf |
|||
Debug.Print "R^2 ="; result(3, 1) |
|||
Debug.Print "F ="; result(4, 1) |
|||
Debug.Print "Degrees of freedom:"; result(4, 2) |
|||
Debug.Print "Standard error of y estimate:"; result(3, 2) |
|||
End Sub</syntaxhighlight>{{out}} |
|||
<pre>coefficients : 1, 2, 3, |
|||
standard errors: 0, 0, 0, |
|||
R^2 = 1 |
|||
F = 7,70461300500498E+31 |
|||
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 |
|||
<syntaxhighlight lang="zkl">var [const] GSL=Import("zklGSL"); // libGSL (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 :=GSL.polyFit(xs,ys,2); |
|||
v.format().println(); |
|||
GSL.Helpers.polyString(v).println(); |
|||
GSL.Helpers.polyEval(v,xs).format().println();</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
1.00,2.00,3.00 |
|||
1 + 2x + 3x^2 |
|||
1.00,6.00,17.00,34.00,57.00,86.00,121.00,162.00,209.00,262.00,321.00 |
|||
</pre> |
|||
Or, using lists: |
|||
{{trans|Common Lisp}} |
|||
Uses the code from [[Multiple regression#zkl]]. |
|||
Example: |
|||
<syntaxhighlight lang="zkl">polyfit(T(T(0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0)), |
|||
T(T(1.0,6.0,17.0,34.0,57.0,86.0,121.0,162.0,209.0,262.0,321.0)), 2) |
|||
.flatten().println();</syntaxhighlight> |
|||
{{out}}<pre>L(1,2,3)</pre> |
Latest revision as of 05:01, 27 April 2024
Find an approximating polynomial of known degree for a given data.
You are encouraged to solve this task according to the task description, using any language you may know.
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 polynomial is:
3 x2 + 2 x + 1
Here, the polynomial's coefficients are (3, 2, 1).
This task is intended as a subtask for Measure relative performance of sorting algorithms implementations.
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)
- Output:
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
Ada
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
function Fit (X, Y : Real_Vector; N : Positive) return Real_Vector is
A : Real_Matrix (0..N, X'Range); -- The plane
begin
for I in A'Range (2) loop
for J in A'Range (1) loop
A (J, I) := X (I)**J;
end loop;
end loop;
return Solve (A * Transpose (A), A * Y);
end Fit;
The function Fit implements least squares approximation of a function defined in the points as specified by the arrays xi and yi. The basis φj is xj, j=0,1,..,N. The implementation is straightforward. First the plane matrix A is created. 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
with Fit;
with Ada.Float_Text_IO; use Ada.Float_Text_IO;
procedure Fitting is
C : constant Real_Vector :=
Fit
( (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0),
(1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0),
2
);
begin
Put (C (0), Aft => 3, Exp => 0);
Put (C (1), Aft => 3, Exp => 0);
Put (C (2), Aft => 3, Exp => 0);
end Fitting;
- Output:
1.000 2.000 3.000
ALGOL 68
MODE FIELD = REAL;
MODE
VEC = [0]FIELD,
MAT = [0,0]FIELD;
PROC VOID raise index error := VOID: (
print(("stop", new line));
stop
);
COMMENT from http://rosettacode.org/wiki/Matrix_Transpose#ALGOL_68 END COMMENT
OP ZIP = ([,]FIELD in)[,]FIELD:(
[2 LWB in:2 UPB in,1 LWB in:1UPB in]FIELD out;
FOR i FROM LWB in TO UPB in DO
out[,i]:=in[i,]
OD;
out
);
COMMENT from http://rosettacode.org/wiki/Matrix_multiplication#ALGOL_68 END COMMENT
OP * = (VEC a,b)FIELD: ( # basically the dot product #
FIELD result:=0;
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI;
FOR i FROM LWB a TO UPB a DO result+:= a[i]*b[i] OD;
result
);
OP * = (VEC a, MAT b)VEC: ( # overload vector times matrix #
[2 LWB b:2 UPB b]FIELD result;
IF LWB a/=LWB b OR UPB a/=UPB b THEN raise index error FI;
FOR j FROM 2 LWB b TO 2 UPB b DO result[j]:=a*b[,j] OD;
result
);
OP * = (MAT a, b)MAT: ( # overload matrix times matrix #
[LWB a:UPB a, 2 LWB b:2 UPB b]FIELD result;
IF 2 LWB a/=LWB b OR 2 UPB a/=UPB b THEN raise index error FI;
FOR k FROM LWB result TO UPB result DO result[k,]:=a[k,]*b OD;
result
);
COMMENT from http://rosettacode.org/wiki/Pyramid_of_numbers#ALGOL_68 END COMMENT
OP / = (VEC a, MAT b)VEC: ( # vector division #
[LWB a:UPB a,1]FIELD transpose a;
transpose a[,1]:=a;
(transpose a/b)[,1]
);
OP / = (MAT a, MAT b)MAT:( # matrix division #
[LWB b:UPB b]INT p ;
INT sign;
[,]FIELD lu = lu decomp(b, p, sign);
[LWB a:UPB a, 2 LWB a:2 UPB a]FIELD out;
FOR col FROM 2 LWB a TO 2 UPB a DO
out[,col] := lu solve(b, lu, p, a[,col]) [@LWB out[,col]]
OD;
out
);
FORMAT int repr = $g(0)$,
real repr = $g(-7,4)$;
PROC fit = (VEC x, y, INT order)VEC:
BEGIN
[0:order, LWB x:UPB x]FIELD a; # the plane #
FOR i FROM 2 LWB a TO 2 UPB a DO
FOR j FROM LWB a TO UPB a DO
a [j, i] := x [i]**j
OD
OD;
( y * ZIP a ) / ( a * ZIP a )
END # fit #;
PROC print polynomial = (VEC x)VOID: (
BOOL empty := TRUE;
FOR i FROM UPB x BY -1 TO LWB x DO
IF x[i] NE 0 THEN
IF x[i] > 0 AND NOT empty THEN print ("+") FI;
empty := FALSE;
IF x[i] NE 1 OR i=0 THEN
IF ENTIER x[i] = x[i] THEN
printf((int repr, x[i]))
ELSE
printf((real repr, x[i]))
FI
FI;
CASE i+1 IN
SKIP,print(("x"))
OUT
printf(($"x**"g(0)$,i))
ESAC
FI
OD;
IF empty THEN print("0") FI;
print(new line)
);
fitting: BEGIN
VEC c =
fit
( (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0),
(1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0),
2
);
print polynomial(c);
VEC d =
fit
( (0, 1, 2, 3, 4, 5, 6, 7, 8, 9),
(2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0),
2
);
print polynomial(d)
END # fitting #
- Output:
3x**2+2x+1 1.0848x**2+10.3552x-0.6164
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
}
Examples:
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
- Output:
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
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;
}
- Output:
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
BBC BASIC
The code listed below is good for up to 10000 data points and fits an order-5 polynomial, so the test data for this task is hardly challenging!
INSTALL @lib$+"ARRAYLIB"
Max% = 10000
DIM vector(5), matrix(5,5)
DIM x(Max%), x2(Max%), x3(Max%), x4(Max%), x5(Max%)
DIM x6(Max%), x7(Max%), x8(Max%), x9(Max%), x10(Max%)
DIM y(Max%), xy(Max%), x2y(Max%), x3y(Max%), x4y(Max%), x5y(Max%)
npts% = 11
x() = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
y() = 1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321
sum_x = SUM(x())
x2() = x() * x() : sum_x2 = SUM(x2())
x3() = x() * x2() : sum_x3 = SUM(x3())
x4() = x2() * x2() : sum_x4 = SUM(x4())
x5() = x2() * x3() : sum_x5 = SUM(x5())
x6() = x3() * x3() : sum_x6 = SUM(x6())
x7() = x3() * x4() : sum_x7 = SUM(x7())
x8() = x4() * x4() : sum_x8 = SUM(x8())
x9() = x4() * x5() : sum_x9 = SUM(x9())
x10() = x5() * x5() : sum_x10 = SUM(x10())
sum_y = SUM(y())
xy() = x() * y() : sum_xy = SUM(xy())
x2y() = x2() * y() : sum_x2y = SUM(x2y())
x3y() = x3() * y() : sum_x3y = SUM(x3y())
x4y() = x4() * y() : sum_x4y = SUM(x4y())
x5y() = x5() * y() : sum_x5y = SUM(x5y())
matrix() = \
\ npts%, sum_x, sum_x2, sum_x3, sum_x4, sum_x5, \
\ sum_x, sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, \
\ sum_x2, sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, \
\ sum_x3, sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, \
\ sum_x4, sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, \
\ sum_x5, sum_x6, sum_x7, sum_x8, sum_x9, sum_x10
vector() = \
\ sum_y, sum_xy, sum_x2y, sum_x3y, sum_x4y, sum_x5y
PROC_invert(matrix())
vector() = matrix().vector()
@% = &2040A
PRINT "Polynomial coefficients = "
FOR term% = 5 TO 0 STEP -1
PRINT ;vector(term%) " * x^" STR$(term%)
NEXT
- Output:
Polynomial coefficients = 0.0000 * x^5 -0.0000 * x^4 0.0002 * x^3 2.9993 * x^2 2.0012 * x^1 0.9998 * x^0
C
Include file (to make the code reusable easily) named polifitgsl.h
#ifndef _POLIFITGSL_H
#define _POLIFITGSL_H
#include <gsl/gsl_multifit.h>
#include <stdbool.h>
#include <math.h>
bool polynomialfit(int obs, int degree,
double *dx, double *dy, double *store); /* n, p */
#endif
Implementation (the examples here helped alot to code this quickly):
#include "polifitgsl.h"
bool polynomialfit(int obs, int degree,
double *dx, double *dy, double *store) /* n, p */
{
gsl_multifit_linear_workspace *ws;
gsl_matrix *cov, *X;
gsl_vector *y, *c;
double chisq;
int i, j;
X = gsl_matrix_alloc(obs, degree);
y = gsl_vector_alloc(obs);
c = gsl_vector_alloc(degree);
cov = gsl_matrix_alloc(degree, degree);
for(i=0; i < obs; i++) {
for(j=0; j < degree; j++) {
gsl_matrix_set(X, i, j, pow(dx[i], j));
}
gsl_vector_set(y, i, dy[i]);
}
ws = gsl_multifit_linear_alloc(obs, degree);
gsl_multifit_linear(X, y, c, cov, &chisq, ws);
/* store result ... */
for(i=0; i < degree; i++)
{
store[i] = gsl_vector_get(c, i);
}
gsl_multifit_linear_free(ws);
gsl_matrix_free(X);
gsl_matrix_free(cov);
gsl_vector_free(y);
gsl_vector_free(c);
return true; /* we do not "analyse" the result (cov matrix mainly)
to know if the fit is "good" */
}
Testing:
#include <stdio.h>
#include "polifitgsl.h"
#define NP 11
double x[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
double y[] = {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321};
#define DEGREE 3
double coeff[DEGREE];
int main()
{
int i;
polynomialfit(NP, DEGREE, x, y, coeff);
for(i=0; i < DEGREE; i++) {
printf("%lf\n", coeff[i]);
}
return 0;
}
- Output:
1.000000 2.000000 3.000000
C#
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<double> qr = v.QR();
// Math.Net doesn't have an "economy" QR, so:
// cut R short to square upper triangle, then recompute Q
var r = qr.R.SubMatrix(0, degree + 1, 0, degree + 1);
var q = v.Multiply(r.Inverse());
var p = r.Inverse().Multiply(q.TransposeThisAndMultiply(yv));
return p.Column(0).ToArray();
}
Example:
static void Main(string[] args)
{
const int degree = 2;
var x = new[] {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
var y = new[] {1.0, 6.0, 17.0, 34.0, 57.0, 86.0, 121.0, 162.0, 209.0, 262.0, 321.0};
var p = Polyfit(x, y, degree);
foreach (var d in p) Console.Write("{0} ",d);
Console.WriteLine();
for (int i = 0; i < x.Length; i++ )
Console.WriteLine("{0} => {1} diff {2}", x[i], Polynomial.Evaluate(x[i], p), y[i] - Polynomial.Evaluate(x[i], p));
Console.ReadKey(true);
}
C++
#include <algorithm>
#include <iostream>
#include <numeric>
#include <vector>
void polyRegression(const std::vector<int>& x, const std::vector<int>& y) {
int n = x.size();
std::vector<int> r(n);
std::iota(r.begin(), r.end(), 0);
double xm = std::accumulate(x.begin(), x.end(), 0.0) / x.size();
double ym = std::accumulate(y.begin(), y.end(), 0.0) / y.size();
double x2m = std::transform_reduce(r.begin(), r.end(), 0.0, std::plus<double>{}, [](double a) {return a * a; }) / r.size();
double x3m = std::transform_reduce(r.begin(), r.end(), 0.0, std::plus<double>{}, [](double a) {return a * a * a; }) / r.size();
double x4m = std::transform_reduce(r.begin(), r.end(), 0.0, std::plus<double>{}, [](double a) {return a * a * a * a; }) / r.size();
double xym = std::transform_reduce(x.begin(), x.end(), y.begin(), 0.0, std::plus<double>{}, std::multiplies<double>{});
xym /= fmin(x.size(), y.size());
double x2ym = std::transform_reduce(x.begin(), x.end(), y.begin(), 0.0, std::plus<double>{}, [](double a, double b) { return a * a * b; });
x2ym /= fmin(x.size(), y.size());
double sxx = x2m - xm * xm;
double sxy = xym - xm * ym;
double sxx2 = x3m - xm * x2m;
double sx2x2 = x4m - x2m * x2m;
double sx2y = x2ym - x2m * ym;
double b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
double c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
double a = ym - b * xm - c * x2m;
auto abc = [a, b, c](int xx) {
return a + b * xx + c * xx*xx;
};
std::cout << "y = " << a << " + " << b << "x + " << c << "x^2" << std::endl;
std::cout << " Input Approximation" << std::endl;
std::cout << " x y y1" << std::endl;
auto xit = x.cbegin();
auto xend = x.cend();
auto yit = y.cbegin();
auto yend = y.cend();
while (xit != xend && yit != yend) {
printf("%2d %3d %5.1f\n", *xit, *yit, abc(*xit));
xit = std::next(xit);
yit = std::next(yit);
}
}
int main() {
using namespace std;
vector<int> x(11);
iota(x.begin(), x.end(), 0);
vector<int> y{ 1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321 };
polyRegression(x, y);
return 0;
}
- Output:
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
Common Lisp
Uses the routine (lsqr A b) from Multiple regression and (mtp A) from Matrix transposition.
;; Least square fit of a polynomial of order n the x-y-curve.
(defun polyfit (x y n)
(let* ((m (cadr (array-dimensions x)))
(A (make-array `(,m ,(+ n 1)) :initial-element 0)))
(loop for i from 0 to (- m 1) do
(loop for j from 0 to n do
(setf (aref A i j)
(expt (aref x 0 i) j))))
(lsqr A (mtp y))))
Example:
(let ((x (make-array '(1 11) :initial-contents '((0 1 2 3 4 5 6 7 8 9 10))))
(y (make-array '(1 11) :initial-contents '((1 6 17 34 57 86 121 162 209 262 321)))))
(polyfit x y 2))
#2A((0.9999999999999759d0) (2.000000000000005d0) (3.0d0))
D
import std.algorithm;
import std.range;
import std.stdio;
auto average(R)(R r) {
auto t = r.fold!("a+b", "a+1")(0, 0);
return cast(double) t[0] / t[1];
}
void polyRegression(int[] x, int[] y) {
auto n = x.length;
auto r = iota(0, n).array;
auto xm = x.average();
auto ym = y.average();
auto x2m = r.map!"a*a".average();
auto x3m = r.map!"a*a*a".average();
auto x4m = r.map!"a*a*a*a".average();
auto xym = x.zip(y).map!"a[0]*a[1]".average();
auto x2ym = x.zip(y).map!"a[0]*a[0]*a[1]".average();
auto sxx = x2m - xm * xm;
auto sxy = xym - xm * ym;
auto sxx2 = x3m - xm * x2m;
auto sx2x2 = x4m - x2m * x2m;
auto sx2y = x2ym - x2m * ym;
auto b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
auto c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
auto a = ym - b * xm - c * x2m;
real abc(int xx) {
return a + b * xx + c * xx * xx;
}
writeln("y = ", a, " + ", b, "x + ", c, "x^2");
writeln(" Input Approximation");
writeln(" x y y1");
foreach (i; 0..n) {
writefln("%2d %3d %5.1f", x[i], y[i], abc(x[i]));
}
}
void main() {
auto x = iota(0, 11).array;
auto y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
polyRegression(x, y);
}
- Output:
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
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[]
Emacs 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)))
- Output:
"3. x^2 + 1.99999999996 x + 1.00000000006"
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
Example
program PolynomalFitting
use fitting
implicit none
! let us test it
integer, parameter :: degree = 2
integer, parameter :: dp = selected_real_kind(15, 307)
integer :: i
real(dp), dimension(11) :: x = (/ (i,i=0,10) /)
real(dp), dimension(11) :: y = (/ 1, 6, 17, 34, &
57, 86, 121, 162, &
209, 262, 321 /)
real(dp), dimension(degree+1) :: a
a = polyfit(x, y, degree)
write (*, '(F9.4)') a
end program
- Output:
(lower powers first, so this seems the opposite of the Python output)
1.0000 2.0000 3.0000
FreeBASIC
General regressions for different polynomials, here it is for degree 2, (3 terms).
#Include "crt.bi" 'for rounding only
Type vector
Dim As Double element(Any)
End Type
Type matrix
Dim As Double element(Any,Any)
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
Next z
Function = ans
End Function
Function matrix.inverse() As matrix
Var ub1=Ubound(this.element,1),ub2=Ubound(this.element,2)
Dim As matrix ans
Dim As vector temp,null_
Redim temp.element(1 To ub1):Redim null_.element(1 To ub1)
Redim ans.element(1 To ub1,1 To ub2)
For a As Integer=1 To ub1
temp=null_
temp.element(a)=1
temp=GaussJordan(temp)
For b As Integer=1 To ub1
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
- Output:
+1 +2*x +3*x^2
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 ]
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)
Go
Library gonum/matrix
package main
import (
"fmt"
"log"
"gonum.org/v1/gonum/mat"
)
func main() {
var (
x = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
y = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
degree = 2
a = Vandermonde(x, degree+1)
b = mat.NewDense(len(y), 1, y)
c = mat.NewDense(degree+1, 1, nil)
)
var qr mat.QR
qr.Factorize(a)
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, d int) *mat.Dense {
x := mat.NewDense(len(a), d, nil)
for i := range a {
for j, p := 0, 1.0; j < d; j, p = j+1, p*a[i] {
x.Set(i, j, p)
}
}
return x
}
- Output:
⎡1.000⎤ ⎢2.000⎥ ⎣3.000⎦
Library go.matrix
Least squares solution using QR decomposition and package go.matrix.
package main
import (
"fmt"
"github.com/skelterjohn/go.matrix"
)
var xGiven = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
var yGiven = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}
var degree = 2
func main() {
m := len(yGiven)
n := degree + 1
y := matrix.MakeDenseMatrix(yGiven, m, 1)
x := matrix.Zeros(m, n)
for i := 0; i < m; i++ {
ip := float64(1)
for j := 0; j < n; j++ {
x.Set(i, j, ip)
ip *= xGiven[i]
}
}
q, r := x.QR()
qty, err := q.Transpose().Times(y)
if err != nil {
fmt.Println(err)
return
}
c := make([]float64, n)
for i := n - 1; i >= 0; i-- {
c[i] = qty.Get(i, 0)
for j := i + 1; j < n; j++ {
c[i] -= c[j] * r.Get(i, j)
}
c[i] /= r.Get(i, i)
}
fmt.Println(c)
}
- Output:
(lowest order coefficient first)
[0.9999999999999758 2.000000000000015 2.999999999999999]
Haskell
Uses module Matrix.LU from hackageDB DSP
import Data.List
import Data.Array
import Control.Monad
import Control.Arrow
import Matrix.LU
ppoly p x = map (x**) p
polyfit d ry = elems $ solve mat vec where
mat = listArray ((1,1), (d,d)) $ liftM2 concatMap ppoly id [0..fromIntegral $ pred d]
vec = listArray (1,d) $ take d ry
- Output:
in GHCi
*Main> polyfit 3 [1,6,17,34,57,86,121,162,209,262,321]
[1.0,2.0,3.0]
HicEst
REAL :: n=10, x(n), y(n), m=3, p(m)
x = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y = (1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
p = 2 ! initial guess for the polynom's coefficients
SOLVE(NUL=Theory()-y(nr), Unknown=p, DataIdx=nr, Iters=iterations)
WRITE(ClipBoard, Name) p, iterations
FUNCTION Theory()
! called by the solver of the SOLVE function. All variables are global
Theory = p(1)*x(nr)^2 + p(2)*x(nr) + p(3)
END
- Output:
SOLVE performs a (nonlinear) least-square fit (Levenberg-Marquardt): p(1)=2.997135145; p(2)=2.011348347; p(3)=0.9906627242; iterations=19;
Hy
(import [numpy [polyfit]])
(setv x (range 11))
(setv y [1 6 17 34 57 86 121 162 209 262 321])
(print (polyfit x y 2))
J
Y=:1 6 17 34 57 86 121 162 209 262 321
(%. ^/~@x:@i.@#) Y
1 2 3 0 0 0 0 0 0 0 0
Note that this implementation does not use floating point numbers, so we do not introduce floating point errors. Using exact arithmetic has a speed penalty, but for small problems like this it is inconsequential.
The above solution fits a polynomial of order 11 (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:
Y %. (i.3) ^/~ i.#Y
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.)
Java
import java.util.Arrays;
import java.util.function.IntToDoubleFunction;
import java.util.stream.IntStream;
public class PolynomialRegression {
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(x).map(a -> a * a).average().orElse(Double.NaN);
double x3m = Arrays.stream(x).map(a -> a * a * a).average().orElse(Double.NaN);
double x4m = Arrays.stream(x).map(a -> a * a * a * a).average().orElse(Double.NaN);
double xym = 0.0;
for (int i = 0; i < x.length && i < y.length; ++i) {
xym += x[i] * y[i];
}
xym /= Math.min(x.length, y.length);
double x2ym = 0.0;
for (int i = 0; i < x.length && i < y.length; ++i) {
x2ym += x[i] * x[i] * y[i];
}
x2ym /= Math.min(x.length, y.length);
double sxx = x2m - xm * xm;
double sxy = xym - xm * ym;
double sxx2 = x3m - xm * x2m;
double sx2x2 = x4m - x2m * x2m;
double sx2y = x2ym - x2m * ym;
double b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
double c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2);
double a = ym - b * xm - c * x2m;
IntToDoubleFunction abc = (int xx) -> a + b * xx + c * xx * xx;
System.out.println("y = " + a + " + " + b + "x + " + c + "x^2");
System.out.println(" Input Approximation");
System.out.println(" x y y1");
for (int i = 0; i < n; ++i) {
System.out.printf("%2d %3d %5.1f\n", x[i], y[i], abc.applyAsDouble(x[i]));
}
}
public static void main(String[] args) {
int[] x = IntStream.range(0, 11).toArray();
int[] y = new int[]{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321};
polyRegression(x, y);
}
}
- Output:
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
jq
Adapted from 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
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
- Output:
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
Julia
The least-squares fit problem for a degree n can be solved with the built-in backslash operator (coefficients in increasing order of degree):
polyfit(x::Vector, y::Vector, deg::Int) = collect(v ^ p for v in x, p in 0:deg) \ y
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)
- Output:
polyfit(x, y, 2) = [1.0, 2.0, 3.0]
Kotlin
// version 1.1.51
fun polyRegression(x: IntArray, y: IntArray) {
val xm = x.average()
val ym = y.average()
val x2m = x.map { it * it }.average()
val x3m = x.map { it * it * it }.average()
val x4m = x.map { it * it * it * it }.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()
val sxx = x2m - xm * xm
val sxy = xym - xm * ym
val sxx2 = x3m - xm * x2m
val sx2x2 = x4m - x2m * x2m
val sx2y = x2ym - x2m * ym
val b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
val c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
val a = ym - b * xm - c * x2m
fun abc(xx: Int) = a + b * xx + c * xx * xx
println("y = $a + ${b}x + ${c}x^2\n")
println(" Input Approximation")
println(" x y y1")
for ((xi, yi) in x zip y) {
System.out.printf("%2d %3d %5.1f\n", xi, yi, abc(xi))
}
}
fun main() {
val x = IntArray(11) { it }
val y = intArrayOf(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
polyRegression(x, y)
}
- Output:
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
Lua
function eval(a,b,c,x)
return a + (b + c * x) * x
end
function regression(xa,ya)
local n = #xa
local xm = 0.0
local ym = 0.0
local x2m = 0.0
local x3m = 0.0
local x4m = 0.0
local xym = 0.0
local x2ym = 0.0
for i=1,n do
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]
end
xm = xm / n
ym = ym / n
x2m = x2m / n
x3m = x3m / n
x4m = x4m / n
xym = xym / n
x2ym = x2ym / n
local sxx = x2m - xm * xm
local sxy = xym - xm * ym
local sxx2 = x3m - xm * x2m
local sx2x2 = x4m - x2m * x2m
local sx2y = x2ym - x2m * ym
local b = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
local c = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
local a = ym - b * xm - c * x2m
print("y = "..a.." + "..b.."x + "..c.."x^2")
for i=1,n do
print(string.format("%2d %3d %3d", xa[i], ya[i], eval(a, b, c, xa[i])))
end
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}
regression(xa, ya)
- Output:
y = 1 + 2x + 3x^2 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
Maple
with(CurveFitting);
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:
3*x^2+2*x+1
Mathematica/Wolfram Language
Using the built-in "Fit" function.
data = Transpose@{Range[0, 10], {1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}};
Fit[data, {1, x, x^2}, x]
Second version: using built-in "InterpolatingPolynomial" function.
Simplify@InterpolatingPolynomial[{{0, 1}, {1, 6}, {2, 17}, {3, 34}, {4, 57}, {5, 86}, {6, 121}, {7, 162}, {8, 209}, {9, 262}, {10, 321}}, x]
WolframAlpha version:
curve fit (0,1), (1,6), (2,17), (3,34), (4,57), (5,86), (6,121), (7,162), (8,209), (9,262), (10,321)
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.
>> 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
МК-61/52
Part 1:
ПC С/П ПD ИП9 + П9 ИПC ИП5 + П5
ИПC x^2 П2 ИП6 + П6 ИП2 ИПC * ИП7
+ П7 ИП2 x^2 ИП8 + П8 ИПC ИПD *
ИПA + ПA ИП2 ИПD * ИПB + ПB ИПD
КИП4 С/П БП 00
Input: В/О x1 С/П y1 С/П x2 С/П y2 С/П ...
Part 2:
ИП5 ПC ИП6 ПD П2 ИП7 П3 ИП4 ИПD *
ИПC ИП5 * - ПD ИП4 ИП7 * ИПC ИП6
* - П7 ИП4 ИПA * ИПC ИП9 * -
ПA ИП4 ИП3 * ИП2 ИП5 * - П3 ИП4
ИП8 * ИП2 ИП6 * - П8 ИП4 ИПB *
ИП2 ИП9 * - ИПD * ИП3 ИПA * -
ИПD ИП8 * ИП7 ИП3 * - / ПB ИПA
ИПB ИП7 * - ИПD / ПA ИП9 ИПB ИП6
* - ИПA ИП5 * - ИП4 / П9 С/П
Result: Р9 = a0, РA = a1, РB = a2.
Modula-2
MODULE PolynomialRegression;
FROM FormatString IMPORT FormatString;
FROM RealStr IMPORT RealToStr;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;
PROCEDURE Eval(a,b,c,x : REAL) : REAL;
BEGIN
RETURN a + b*x + c*x*x;
END Eval;
PROCEDURE Regression(x,y : ARRAY OF INTEGER);
VAR
n,i : INTEGER;
xm,x2m,x3m,x4m : REAL;
ym : REAL;
xym,x2ym : REAL;
sxx,sxy,sxx2,sx2x2,sx2y : REAL;
a,b,c : REAL;
buf : ARRAY[0..63] OF CHAR;
BEGIN
n := SIZE(x)/SIZE(INTEGER);
xm := 0.0;
ym := 0.0;
x2m := 0.0;
x3m := 0.0;
x4m := 0.0;
xym := 0.0;
x2ym := 0.0;
FOR i:=0 TO n-1 DO
xm := xm + FLOAT(x[i]);
ym := ym + FLOAT(y[i]);
x2m := x2m + FLOAT(x[i]) * FLOAT(x[i]);
x3m := x3m + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]);
x4m := x4m + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(x[i]);
xym := xym + FLOAT(x[i]) * FLOAT(y[i]);
x2ym := x2ym + FLOAT(x[i]) * FLOAT(x[i]) * FLOAT(y[i]);
END;
xm := xm / FLOAT(n);
ym := ym / FLOAT(n);
x2m := x2m / FLOAT(n);
x3m := x3m / FLOAT(n);
x4m := x4m / FLOAT(n);
xym := xym / FLOAT(n);
x2ym := x2ym / FLOAT(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;
WriteString("y = ");
RealToStr(a, buf);
WriteString(buf);
WriteString(" + ");
RealToStr(b, buf);
WriteString(buf);
WriteString("x + ");
RealToStr(c, buf);
WriteString(buf);
WriteString("x^2");
WriteLn;
FOR i:=0 TO n-1 DO
FormatString("%2i %3i ", buf, x[i], y[i]);
WriteString(buf);
RealToStr(Eval(a,b,c,FLOAT(x[i])), buf);
WriteString(buf);
WriteLn;
END;
END Regression;
TYPE R = ARRAY[0..10] OF INTEGER;
VAR
x,y : R;
BEGIN
x := R{0,1,2,3,4,5,6,7,8,9,10};
y := R{1,6,17,34,57,86,121,162,209,262,321};
Regression(x,y);
ReadChar;
END PolynomialRegression.
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)
- Output:
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
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
- Output:
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
Octave
x = [0:10];
y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321];
coeffs = polyfit(x, y, 2)
PARI/GP
Lagrange interpolating polynomial:
polinterpolate([0,1,2,3,4,5,6,7,8,9,10],[1,6,17,34,57,86,121,162,209,262,321])
In newer versions, this can be abbreviated:
polinterpolate([0..10],[1,6,17,34,57,86,121,162,209,262,321])
- Output:
3*x^2 + 2*x + 1
Least-squares fit:
V=[1,6,17,34,57,86,121,162,209,262,321]~;
M=matrix(#V,3,i,j,(i-1)^(j-1));Polrev(matsolve(M~*M,M~*V))
Code thanks to Bill Allombert
- Output:
3*x^2 + 2*x + 1
Least-squares polynomial fit in its own function:
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)
Perl
This code identical to that of Multiple regression task.
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;
- Output:
const 1.000 X 2.000 X**2 3.000
PDL Alternative:
#!/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($_)
}
- Output:
3.00000037788248 * $x**2 + 1.99999750988868 * $x + 1.00000180493936
Phix
You can run this online here.
-- demo\rosetta\Polynomial_regression.exw with javascript_semantics constant x = {0,1,2,3,4,5,6,7,8,9,10}, y = {1,6,17,34,57,86,121,162,209,262,321}, n = length(x) function regression() atom xm = 0, ym = 0, x2m = 0, x3m = 0, x4m = 0, xym = 0, x2ym = 0 for i=1 to n do atom xi = x[i], yi = y[i] xm += xi ym += yi x2m += power(xi,2) x3m += power(xi,3) x4m += power(xi,4) xym += xi*yi x2ym += power(xi,2)*yi end for xm /= n ym /= n x2m /= n x3m /= n x4m /= n xym /= n x2ym /= n atom Sxx = x2m-power(xm,2), Sxy = xym-xm*ym, Sxx2 = x3m-xm*x2m, Sx2x2 = x4m-power(x2m,2), Sx2y = x2ym-x2m*ym, B = (Sxy*Sx2x2-Sx2y*Sxx2)/(Sxx*Sx2x2-power(Sxx2,2)), C = (Sx2y*Sxx-Sxy*Sxx2)/(Sxx*Sx2x2-power(Sxx2,2)), A = ym-B*xm-C*x2m return {C,B,A} end function atom {a,b,c} = regression() function f(atom x) return a*x*x+b*x+c end function printf(1,"y=%gx^2+%gx+%g\n",{a,b,c}) printf(1,"\n x y f(x)\n") for i=1 to n do printf(1," %2d %3d %3g\n",{x[i],y[i],f(x[i])}) end for -- And a simple plot (re-using x,y from above) include pGUI.e include IupGraph.e function get_data(Ihandle graph) integer {w,h} = IupGetIntInt(graph,"DRAWSIZE") IupSetInt(graph,"YTICK",iff(h<240?iff(h<150?80:40):20)) return {{x,y,CD_RED}} end function IupOpen() Ihandle graph = IupGraph(get_data,"RASTERSIZE=640x440") IupSetAttributes(graph,"XTICK=1,XMIN=0,XMAX=10") IupSetAttributes(graph,"YTICK=20,YMIN=0,YMAX=320") Ihandle dlg = IupDialog(graph,`TITLE="simple plot"`) IupSetAttributes(dlg,"MINSIZE=245x150") IupShow(dlg) if platform()!=JS then IupMainLoop() IupClose() end if
- Output:
(plus a simple graphical plot, as per Racket)
y=3x^2+2x+1 x y f(x) 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
PowerShell
function qr([double[][]]$A) {
$m,$n = $A.count, $A[0].count
$pm,$pn = ($m-1), ($n-1)
[double[][]]$Q = 0..($m-1) | foreach{$row = @(0) * $m; $row[$_] = 1; ,$row}
[double[][]]$R = $A | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
foreach ($h in 0..$pn) {
[double[]]$u = $R[$h..$pm] | foreach{$_[$h]}
[double]$nu = $u | foreach {[double]$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)}
$u[0] -= if ($u[0] -lt 1) {$nu} else {-$nu}
[double]$nu = $u | foreach {$sq = 0} {$sq += $_*$_} {[Math]::Sqrt($sq)}
[double[]]$u = $u | foreach { $_/$nu}
[double[][]]$v = 0..($u.Count - 1) | foreach{$i = $_; ,($u | foreach{2*$u[$i]*$_})}
[double[][]]$CR = $R | foreach{$row = $_; ,@(0..$pn | foreach{$row[$_]})}
[double[][]]$CQ = $Q | foreach{$row = $_; ,@(0..$pm | foreach{$row[$_]})}
foreach ($i in $h..$pm) {
foreach ($j in $h..$pn) {
$R[$i][$j] -= $h..$pm | foreach {[double]$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CR[$_][$j]} {$sum}
}
}
if (0 -eq $h) {
foreach ($i in $h..$pm) {
foreach ($j in $h..$pm) {
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i][$_]*$CQ[$_][$j]} {$sum}
}
}
} else {
$p = $h-1
foreach ($i in $h..$pm) {
foreach ($j in 0..$p) {
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
}
foreach ($j in $h..$pm) {
$Q[$i][$j] -= $h..$pm | foreach {$sum = 0} {$sum += $v[$i-$h][$_-$h]*$CQ[$_][$j]} {$sum}
}
}
}
}
foreach ($i in 0..$pm) {
foreach ($j in $i..$pm) {$Q[$i][$j],$Q[$j][$i] = $Q[$j][$i],$Q[$i][$j]}
}
[PSCustomObject]@{"Q" = $Q; "R" = $R}
}
function leastsquares([Double[][]]$A,[Double[]]$y) {
$QR = qr $A
[Double[][]]$Q = $QR.Q
[Double[][]]$R = $QR.R
$m,$n = $A.count, $A[0].count
[Double[]]$z = foreach ($j in 0..($m-1)) {
0..($m-1) | foreach {$sum = 0} {$sum += $Q[$_][$j]*$y[$_]} {$sum}
}
[Double[]]$x = @(0)*$n
for ($i = $n-1; $i -ge 0; $i--) {
for ($j = $i+1; $j -lt $n; $j++) {
$z[$i] -= $x[$j]*$R[$i][$j]
}
$x[$i] = $z[$i]/$R[$i][$i]
}
$x
}
function polyfit([Double[]]$x,[Double[]]$y,$n) {
$m = $x.Count
[Double[][]]$A = 0..($m-1) | foreach{$row = @(1) * ($n+1); ,$row}
for ($i = 0; $i -lt $m; $i++) {
for ($j = $n-1; 0 -le $j; $j--) {
$A[$i][$j] = $A[$i][$j+1]*$x[$i]
}
}
leastsquares $A $y
}
function show($m) {$m | foreach {write-host "$_"}}
$A = @(@(12,-51,4), @(6,167,-68), @(-4,24,-41))
$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^2 X constant"
"$(polyfit $x $y 2)"
- Output:
polyfit X^2 X constant 3 1.99999999999998 1.00000000000005
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.])
Substitute back received coefficients.
>>> yf = numpy.polyval(numpy.poly1d(coeffs), x)
>>> yf
array([ 1., 6., 17., 34., 57., 86., 121., 162., 209., 262., 321.])
Find max absolute error:
>>> '%.1g' % max(y-yf)
'1e-013'
Example
For input arrays `x' and `y':
>>> x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> y = [2.7, 2.8, 31.4, 38.1, 58.0, 76.2, 100.5, 130.0, 149.3, 180.0]
>>> p = numpy.poly1d(numpy.polyfit(x, y, deg=2), variable='N')
>>> print p
2
1.085 N + 10.36 N - 0.6164
Thus we confirm once more that for already sorted sequences the considered quick sort implementation has quadratic dependence on sequence length (see Example section for Python language on Query Performance page).
R
The easiest (and most robust) approach to solve this in R is to use the base package's lm function which will find the least squares solution via a QR decomposition:
x <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
y <- c(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321)
coef(lm(y ~ x + I(x^2)))
- Output:
(Intercept) x I(x^2) 1 2 3
Alternately, use poly:
coef(lm(y ~ poly(x, 2, raw=T)))
- Output:
(Intercept) poly(x, 2, raw = T)1 poly(x, 2, raw = T)2 1 2 3
Racket
#lang racket
(require math plot)
(define xs '(0 1 2 3 4 5 6 7 8 9 10))
(define ys '(1 6 17 34 57 86 121 162 209 262 321))
(define (fit x y n)
(define Y (->col-matrix y))
(define V (vandermonde-matrix x (+ n 1)))
(define VT (matrix-transpose V))
(matrix->vector (matrix-solve (matrix* VT V) (matrix* VT Y))))
(define ((poly v) x)
(for/sum ([c v] [i (in-naturals)])
(* c (expt x i))))
(plot (list (points (map vector xs ys))
(function (poly (fit xs ys 2)))))
- Output:
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 : , où on cherche . On considère et on répartit chaque équation sur chaque dimension:
Posons alors :
Le système d'équations devient : .
D'où :
use MultiVector;
constant @x1 = <0 1 2 3 4 5 6 7 8 9 10>;
constant @y = <1 6 17 34 57 86 121 162 209 262 321>;
constant $x0 = [+] @e[^@x1];
constant $x1 = [+] @x1 Z* @e;
constant $x2 = [+] @x1 »**» 2 Z* @e;
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);
- Output:
1 2 3
REXX
/* REXX ---------------------------------------------------------------
* Implementation of http://keisan.casio.com/exec/system/14059932254941
*--------------------------------------------------------------------*/
xl='0 1 2 3 4 5 6 7 8 9 10'
yl='1 6 17 34 57 86 121 162 209 262 321'
n=11
Do i=1 To n
Parse Var xl x.i xl
Parse Var yl y.i yl
End
xm=0
ym=0
x2m=0
x3m=0
x4m=0
xym=0
x2ym=0
Do i=1 To n
xm=xm+x.i
ym=ym+y.i
x2m=x2m+x.i**2
x3m=x3m+x.i**3
x4m=x4m+x.i**4
xym=xym+x.i*y.i
x2ym=x2ym+(x.i**2)*y.i
End
xm =xm /n
ym =ym /n
x2m=x2m/n
x3m=x3m/n
x4m=x4m/n
xym=xym/n
x2ym=x2ym/n
Sxx=x2m-xm**2
Sxy=xym-xm*ym
Sxx2=x3m-xm*x2m
Sx2x2=x4m-x2m**2
Sx2y=x2ym-x2m*ym
B=(Sxy*Sx2x2-Sx2y*Sxx2)/(Sxx*Sx2x2-Sxx2**2)
C=(Sx2y*Sxx-Sxy*Sxx2)/(Sxx*Sx2x2-Sxx2**2)
A=ym-B*xm-C*x2m
Say 'y='a'+'||b'*x+'c'*x**2'
Say ' Input "Approximation"'
Say ' x y y1'
Do i=1 To 11
Say right(x.i,2) right(y.i,3) format(fun(x.i),5,3)
End
Exit
fun:
Parse Arg x
Return a+b*x+c*x**2
- Output:
y=1+2*x+3*x**2 Input "Approximation" x y y1 0 1 1.000 1 6 6.000 2 17 17.000 3 34 34.000 4 57 57.000 5 86 86.000 6 121 121.000 7 162 162.000 8 209 209.000 9 262 262.000 10 321 321.000
RPL
≪ 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 * / @ the following lines convert the resulting vector into a polynomial equation DUP 'x' STO 1 GET 2 x SIZE FOR j 'X' * x j GET + NEXT EXPAN COLCT ≫ ≫ 'FIT' STO
[1 2 3 4 5 6 7 8 9 10] [1 6 17 34 57 86 121 162 209 262 321] 2 FIT
- Output:
1: '3+2*X+1*X^2'
Ruby
require 'matrix'
def regress x, y, degree
x_data = x.map { |xi| (0..degree).map { |pow| (xi**pow).to_r } }
mx = Matrix[*x_data]
my = Matrix.column_vector(y)
((mx.t * mx).inv * mx.t * my).transpose.to_a[0].map(&:to_f)
end
Testing:
p regress([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
[1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321],
2)
- Output:
[1.0, 2.0, 3.0]
Scala
- Output:
See it yourself by running in your browser Scastie (remote JVM).
object PolynomialRegression extends App {
private def xy = Seq(1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321).zipWithIndex.map(_.swap)
private def polyRegression(xy: Seq[(Int, Int)]): Unit = {
val r = xy.indices
def average[U](ts: Iterable[U])(implicit num: Numeric[U]) = num.toDouble(ts.sum) / ts.size
def x3m: Double = average(r.map(a => a * a * a))
def x4m: Double = average(r.map(a => a * a * a * a))
def x2ym = xy.reduce((a, x) => (a._1 + x._1 * x._1 * x._2, 0))._1.toDouble / xy.size
def xym = xy.reduce((a, x) => (a._1 + x._1 * x._2, 0))._1.toDouble / xy.size
val x2m: Double = average(r.map(a => a * a))
val (xm, ym) = (average(xy.map(_._1)), average(xy.map(_._2)))
val (sxx, sxy) = (x2m - xm * xm, xym - xm * ym)
val sxx2: Double = x3m - xm * x2m
val sx2x2: Double = x4m - x2m * x2m
val sx2y: Double = x2ym - x2m * ym
val c: Double = (sx2y * sxx - sxy * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
val b: Double = (sxy * sx2x2 - sx2y * sxx2) / (sxx * sx2x2 - sxx2 * sxx2)
val a: Double = ym - b * xm - c * x2m
def abc(xx: Int) = a + b * xx + c * xx * xx
println(s"y = $a + ${b}x + ${c}x^2")
println(" Input Approximation")
println(" x y y1")
xy.foreach {el => println(f"${el._1}%2d ${el._2}%3d ${abc(el._1)}%5.1f")}
}
polyRegression(xy)
}
Sidef
func regress(x, y, degree) {
var A = Matrix.build(x.len, degree+1, {|i,j|
x[i]**j
})
var B = Matrix.column_vector(y...)
((A.transpose * A)**(-1) * A.transpose * B).transpose[0]
}
func poly(x) {
3*x**2 + 2*x + 1
}
var coeff = regress(
10.of { _ },
10.of { poly(_) },
2
)
say coeff
- Output:
[1, 2, 3]
Stata
See Factor variables in Stata help for explanations on the c.x##c.x syntax.
. clear
. input x y
0 1
1 6
2 17
3 34
4 57
5 86
6 121
7 162
8 209
9 262
10 321
end
. regress y c.x##c.x
Source | SS df MS Number of obs = 11
-------------+---------------------------------- F(2, 8) = .
Model | 120362 2 60181 Prob > F = .
Residual | 0 8 0 R-squared = 1.0000
-------------+---------------------------------- Adj R-squared = 1.0000
Total | 120362 10 12036.2 Root MSE = 0
------------------------------------------------------------------------------
y | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
x | 2 . . . . .
|
c.x#c.x | 3 . . . . .
|
_cons | 1 . . . . .
------------------------------------------------------------------------------
Swift
let x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
let y = [1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321]
func average(_ input: [Int]) -> Int {
return input.reduce(0, +) / input.count
}
func polyRegression(x: [Int], y: [Int]) {
let xm = average(x)
let ym = average(y)
let x2m = average(x.map { $0 * $0 })
let x3m = average(x.map { $0 * $0 * $0 })
let x4m = average(x.map { $0 * $0 * $0 * $0 })
let xym = average(zip(x,y).map { $0 * $1 })
let x2ym = average(zip(x,y).map { $0 * $0 * $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(xx: Int) -> Int {
return (a + b * xx) + (c * xx * xx)
}
print("y = \(a) + \(b)x + \(c)x^2\n")
print(" Input Approximation")
print(" x y y1")
for i in 0 ..< x.count {
let result = Double(abc(xx: i))
print(String(format: "%2d %3d %5.1f", x[i], y[i], result))
}
}
polyRegression(x: x, y: y)
- Output:
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
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
This will print:
1.0000000000000207 1.9999999999999958 3.0
which is a close approximation to the correct solution.
TI-83 BASIC
DelVar X
seq(X,X,0,10) → L1
{1,6,17,34,57,86,121,162,209,262,321} → L2
QuadReg L1,L2
- Output:
y=ax²+bx+c a=3 b=2 c=1
TI-89 BASIC
DelVar x
seq(x,x,0,10) → xs
{1,6,17,34,57,86,121,162,209,262,321} → ys
QuadReg xs,ys
Disp regeq(x)
seq(expr,var,low,high)
evaluates expr with var bound to integers from low to high and returns a list of the results. →
is the assignment operator.
QuadReg
, "quadratic regression", does the fit and stores the details in a number of standard variables, including regeq, which receives the fitted quadratic (polynomial) function itself.
We then apply that function to the (undefined as ensured by DelVar
) variable x to obtain the expression in terms of x, and display it.
- Output:
3.·x2 + 2.·x + 1.
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.
#import std
#import nat
#import flo
(fit "n") ("x","y") = ..dgelsd\"y" (gang \/*pow float*x iota successor "n")* "x"
test program:
x = <0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.>
y = <1.,6.,17.,34.,57.,86.,121.,162.,209.,262.,321.>
#cast %eL
example = fit2(x,y)
- Output:
<3.000000e+00,2.000000e+00,1.000000e+00>
VBA
Excel VBA has built in capability for line estimation.
Option Base 1
Private Function polynomial_regression(y As Variant, x As Variant, degree As Integer) As Variant
Dim a() As Double
ReDim a(UBound(x), 2)
For i = 1 To UBound(x)
For j = 1 To degree
a(i, j) = x(i) ^ j
Next j
Next i
polynomial_regression = WorksheetFunction.LinEst(WorksheetFunction.Transpose(y), a, True, True)
End Function
Public Sub main()
x = [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}]
y = [{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321}]
result = polynomial_regression(y, x, 2)
Debug.Print "coefficients : ";
For i = UBound(result, 2) To 1 Step -1
Debug.Print Format(result(1, i), "0.#####"),
Next i
Debug.Print
Debug.Print "standard errors: ";
For i = UBound(result, 2) To 1 Step -1
Debug.Print Format(result(2, i), "0.#####"),
Next i
Debug.Print vbCrLf
Debug.Print "R^2 ="; result(3, 1)
Debug.Print "F ="; result(4, 1)
Debug.Print "Degrees of freedom:"; result(4, 2)
Debug.Print "Standard error of y estimate:"; result(3, 2)
End Sub
- Output:
coefficients : 1, 2, 3, standard errors: 0, 0, 0, R^2 = 1 F = 7,70461300500498E+31 Degrees of freedom: 8 Standard error of y estimate: 2,79482284961344E-14
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)
- Output:
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
zkl
Using the GNU Scientific Library
var [const] GSL=Import("zklGSL"); // libGSL (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 :=GSL.polyFit(xs,ys,2);
v.format().println();
GSL.Helpers.polyString(v).println();
GSL.Helpers.polyEval(v,xs).format().println();
- Output:
1.00,2.00,3.00 1 + 2x + 3x^2 1.00,6.00,17.00,34.00,57.00,86.00,121.00,162.00,209.00,262.00,321.00
Or, using lists:
Uses the code from Multiple regression#zkl.
Example:
polyfit(T(T(0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0)),
T(T(1.0,6.0,17.0,34.0,57.0,86.0,121.0,162.0,209.0,262.0,321.0)), 2)
.flatten().println();
- Output:
L(1,2,3)