МНК алгоритм

Метод наименьших квадратов

Задача написать алгоритм для нахождения коэффициентов полинома N степени с помощью МНК по M точек. То есть даны 10 точек, и сказано аппроксимировать для полинома 3-й степени:

y=ax^3 + a1x^2 + a2x+a3

Для этого полинома надо найти a, a1, a2, a3; если для 1 и 2 степени можно написать вручную, то даже для 3-й степени придётся решать систему 4-х уравнений с 4 неизвестными.

Вот тут можно посмотреть, как выглядит аппроксимация для разных функций с помощью мнк.

Как запрограммировать алгоритм, который будет составлять СЛОУ и решать её?

Язык c++.

Обновление

Как в составлении так и в решении. Начал решать пока что такая система получается для полинома N степени. Ex сумма от 1 до M (кол-во точек).

a1*Ex^n     + a2*Ex^n - 1  + ... + an*Ex^0=Eyx^0
a1*Ex^n + 1 + a2*Ex^n      + ... + an*Ex^1=Eyx^1

a1*Ex^2n    + a2*Ex^2n - 1 + ... + an*Ex^n=Eyx^n

Даже после того как перепишу все известные Ex,Ey,Exy... в массив, то я незнаю как решать такую СЛАУ

Решение

Прошу меня простить за кривые отступы. Функция которая находит коэф. полинома N степени

    vector<double>* mnk(vector<vector<double>> pts,int n,int n2)  
    {
        //матрица хранящая неизвестные при нужных нам коэффициентах 
        vector<vector<double>> linSys(n + 1);
double tmp=0;
vector<double> subsid(3*n + 2);//почему 3n + 2? да потому что http://mathhelpplanet.com/static.php?p=onlayn - mnk - i - regressionniy - analiz
for (int i=0;i<2*n + 1;i ++ )//этот цикл формирует первую половину (точнее 2/3) списка величин x^0 x^1...x^2n + 1
{
    for (int j=0;j<n2;j ++ )
    {
        tmp + =pow(pts[j][0],i);
    }
    subsid[i]=tmp;
    tmp=0;
}
for (int i=2*n + 1;i<3*n + 2;i ++ )//этот цикл формирует вторую половину  списка величин x^0*y x^1*y...x^n*y
{
    for (int j=0;j<n2;j ++ )
    {
        tmp + =pow(pts[j][0],i - 2*n - 1)*pts[j][1];
    }
    subsid[i]=tmp;
    tmp=0;
}
//тут уже идёт формирование СЛАУ получаем матрицу linSys
//a1*Ex^n      +  a2*Ex^n  -  1   +  ...  +  an*Ex^0=Eyx^0
//a1*Ex^n  +  1  +  a2*Ex^n       +  ...  +  an*Ex^1=Eyx^1
//  ...            ...         ...       ... 
//a1*Ex^2n     +  a2*Ex^2n  -  1  +  ...  +  an*Ex^n=Eyx^n
//Будем решать СЛАУ в виде A*x=b
//в моем случае linSys*x=b
vector<double>*b=new vector<double>;
b - >resize(n + 1);
for(int i = 0; i < n  +  1; i ++ )
{
    linSys[i].resize(n + 1);
    (*b)[i]=subsid[2*n + 1 + i];
    for(int j=0;j<n + 1;j ++ )
    {
        linSys[i][j]=subsid[n + i - j];
    }
}
vector<double>*x=new vector<double>;
x - >resize(n + 1);
//это решает СЛАУ и пишет результ в x
Gauss(linSys,*b,x,n + 1);
return x;

}

Функция которая решает СЛАУ

#include <math.h>
#include <iostream>
#include <vector>
#include <stack>
//результат записывает в x, СЛАУ представима как A*x=B
void Gauss(vector<vector<double>> A, vector<double> b, vector<double> *x, int n)
    {
std::stack<std::pair<size_t,size_t> > swaps;
int i, j, k, t;
double kof, s;
for (i = n  -  1; i > 0;  -- i) {
    for (t=i, j=i - 1; j >= 0;  -- j) {//Поиск максимума в столбце
        if (fabs(A[i][t]) < fabs(A[i][j])) {
            t = j;
        }
    }
    if (A[i][t] == 0.0) {
        return;
    }
    if (t != i) {
        for (k = n  -  1; k >= 0;  -- k) {
            std::swap(A[k][t],A[k][i]);
        }
        swaps.push(std::pair<size_t, size_t>(i,t));
    }
    for (j=i - 1; j>=0;  -- j) {//
        kof = A[j][i] / A[i][i];//
        A[j][i] = 0.0;//
        b[j]  - = b[i] * kof;
        for (k = i  -  1; k >= 0;  -- k) {
            A[j][k]  - = A[i][k] * kof;
        }
    }
}
for (i = 0; i < n;  ++ i)
{
    s = 0.0;
    for (j = i  -  1; j >= 0;  -- j) {
        s  + = A[i][j] * (*x)[j];
    }
    (*x)[i] = (b[i]  -  s) / A[i][i];
}
while (!swaps.empty()) {
    //int j=swaps.top().first;
//  int y=swaps.top().second;
    std::swap((*x)[swaps.top().first], (*x)[swaps.top().second]);
    swaps.pop();
}

}