Polynomial regression: Difference between revisions

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

{{trans|Java}}

<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;

}</lang>

{{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|C#|C sharp}}==

{{libheader|Math.Net}}

=={{header|Common Lisp}}==

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

<lang lisp>;; Least square fit of a polynomial of order n the x-y-curve.

(defun polyfit (x y n)

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

(A (make-array `(,m ,(+ n 1)) :initial-element 0)))

(loop for i from 0 to (- m 1) do

(loop for j from 0 to n do

(setf (aref A i j)

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

(lsqr A (mtp y))))</lang>

Example:

<lang lisp>(let ((x (make-array '(1 11) :initial-contents '((0 1 2 3 4 5 6 7 8 9 10))))

(y (make-array '(1 11) :initial-contents '((1 6 17 34 57 86 121 162 209 262 321)))))

(polyfit x y 2))

#2A((0.9999999999999759d0) (2.000000000000005d0) (3.0d0))</lang>