Аппроксимация степенным полиномом используя метод наименьших квадратов

// approximation.cpp: определяет точку входа для консольного приложения.
#include <fstream>
#include <math.h>
#include <conio.h>
#include <windows.h>
#include <stdio.h>
#include <iostream>
using namespace std;
double tabl_X[5] = { 0, 20, 50, 60, 75 }; //исходные данные X
double tabl_Y[5] = { 0, 2.67, 7.04, 9.30, 10.82 }; //исходные данные Y
const int N = 5; //число точек
const double Emax = 3;//допустимая среднеквадратичная ошибка аппроксимации
int m; //ширина матрицы Х, длина вектора А (m <= N)
int i, j; //счетчики
double** X; //матрица Х - рабочая
double** Xt; //матрица Х - транспонированная
double** XtX; //произведение XtX
double* XtY; //произведение XtY
double** XtXb; //обратная XtX
double* A; //вектор коэффициентов
double e[5]; //вектор квадратичных ошибок
double E; //ошибка аппроксимации
//Функция печати матрицы
void print_matr(double** X, int n, int m, int k) {
cout.precision(k);
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++)
cout << X[i][j] << " ";
cout << endl;
}
}//========================================================
//Функция печати вектора
void print_vect(double* X, int n, int k) {
cout.precision(k); //число выводимых после точки знаков
for (i = 0; i < n; i++)
cout << X[i] << " ";
cout << endl;
}//========================================================
//Функция освобождения памяти из-под матрицы
void free_mem(double** X, int n) {
for (i = 0; i < n; i++)
delete[] X[i];
delete[] X;
}//========================================================
//функция вычеркивания строки и столбца
void get_matr(double** matr, int n, double** temp_matr, int indRow, int indCol) {
int ki = 0, kj;
for (i = 0; i < n; i++) {
if (i != indRow) {
kj = 0;
for (j = 0; j < n; j++) {
if (j != indCol) {
temp_matr[ki][kj] = matr[i][j];
kj++;
}
}
ki++;
}
}
}//========================================================
//функция вычисления определителя матрицы
double determ(double** X, int m) {
double sum = 0;
int n;
if (m == 1)
sum = X[0][0]; //расчет определителя для матрцицы единичной размерноси
else if (m == 2)
sum = X[0][0] * X[1][1] - X[1][0] * X[0][1]; //расчет определителя для матрцицы размерноси 2
else {
for (int i = 0; i < m; i++) { //расчет определителя для матрцицы размерности больше двух
n = m - 1;
double** temp_matr = new double* [n];
for (j = 0; j < n; j++)
temp_matr[j] = new double[n];
get_matr(X, m, temp_matr, 0, i);
double det = determ(temp_matr, n); //Вызывается функция расчета определителя
sum = sum + pow(-1.0, 2 + i) * X[0][i] * det; //расчет непосредственно определителя
free_mem(temp_matr, n); //освобождаем память из под матрицы temp_matr
}
}
return sum;
}//========================================================
int main(int argc, char** argv)
{
m = 1;
do
{
m++; //на первой итерации степень полинома равна 1, соответственно длина вектора А равна 2
//формируем матрицу Х
X = new double* [N];
for (i = 0; i < N; i++)
X[i] = new double[m];
for (i = 0; i < N; i++) {
X[i][0] = 1;
X[i][1] = tabl_X[i];
if (m > 2)
for (j = 2; j < m; j++)
X[i][j] = pow(tabl_X[i], j);
}
//получаем транспонированную матрицу X
Xt = new double* [m];
for (i = 0; i < m; i++)
Xt[i] = new double[N];
for (i = 0; i < m; i++)
for (j = 0; j < N; j++)
Xt[i][j] = X[j][i];
//вычисляем произведения XtX и XtY
XtX = new double* [m];
for (i = 0; i < m; i++)
XtX[i] = new double[m];
XtY = new double[m];
double sum;
int k;
for (i = 0; i < m; i++) {
for (j = 0; j < m; j++) {
sum = 0;
for (k = 0; k < N; k++)
sum = sum + Xt[i][k] * X[k][j]; //перемножение строки Хt на столбец Х
XtX[i][j] = sum;
}
sum = 0;
for (k = 0; k < N; k++)
sum = sum + Xt[i][k] * tabl_Y[k]; //перемножение Xt на вектор Y
XtY[i] = sum;
}
//вычисляем обратную матрицу
XtXb = new double* [m];
for (i = 0; i < m; i++)
XtXb[i] = new double[m];
double det = determ(XtX, m); //Вызывается функция расчета определителя
if (det != 0) {
for (int i = 0; i < m; i++) {
for (int j = 0; j < m; j++) {
int n = m - 1;
double** temp_matr = new double* [n];
for (int k = 0; k < n; k++)
temp_matr[k] = new double[n];
get_matr(XtX, m, temp_matr, i, j);
double det1 = determ(temp_matr, n);
XtXb[i][j] = pow(-1.0, i + j + 2) * det1 / det; //рассчитываем элементы обратной матрицы XtXb
free_mem(temp_matr, n); //освобождаем память из под матрицы temp_matr
}
}
}
else
cout << "Невозможно найти обратную матрицу" << endl;
//вычисляем вектор коэффициентов
A = new double[m];
for (i = 0; i < m; i++) {
sum = 0;
for (k = 0; k < m; k++)
sum = sum + XtXb[i][k] * XtY[k]; //вектор А рассчитывается как произведение матрицы XtXb на вектор XtY
A[i] = sum;
}
//Вычисляем ошибку аппроксимации
for (i = 0; i < N; i++) {
sum = 0;
for (j = 0; j < m; j++)
sum = sum + pow(tabl_X[i], j) * A[j];
e[i] = pow(tabl_Y[i] - sum, 2); //расчет вектора квадратичных ошибок
}
sum = 0;
for (i = 0; i < N; i++)
sum = sum + e[i];
E = pow(sum / (1.0 * N), 0.5); //расчет среднеквадратичной ошибки аппроксимации
} while ((E > Emax) && (m < N));
setlocale(LC_ALL, ".ACP");
printf("Исходные данные:\n\n");
printf(" X | Y \n");
printf("-----------\n");
for (i = 0; i < N; i++)
printf(" %2.0lf | %5.2lf \n", tabl_X[i], tabl_Y[i]);
printf("\nEmax = %.2lf\n\n", Emax);
printf("============================================\n\n");
cout << fixed;
printf("Результаты расчетов:\n\n");
printf("Порядок регрессионного полинома M = %d\n\n", m - 1);
printf("Матрица Х:\n");
print_matr(X, N, m, 0);
printf("\nМатрица Хt - транспонированная Х:\n");
print_matr(Xt, m, N, 0);
printf("\nМатрица ХtX = Xt * X:\n");
print_matr(XtX, m, m, 0);
printf("\nВектор ХtY = Xt * tabl_Y:\n");
print_vect(XtY, m, 3);
printf("\nМатрица ХtYb - матрица, обратная XtY:\n");
print_matr(XtXb, m, m, 7);
printf("\nВектор А - вектор коэффициентов регрессионного полинома:\n");
print_vect(A, m, 9);
printf("\nCреднеквадратичная ошибка аппроксимации E = %.5lf\n", E);
_getch();
return 0;
}
 
 
Листинг программы:
 
// approximation2.cpp: определяет точку входа для консольного приложения.
#include <fstream>
#include <conio.h>
#include <windows.h>
#include <stdio.h>
#include <iostream>
using namespace std;
double tabl_X[5] = { 0, 20, 50, 60, 75 }; //исходные данные X
double tabl_Y[5] = { 0, 2.67, 7.04, 9.30, 10.82 }; //исходные данные Y
const int N = 5; //число точек
int m; //длина вектора А (m <= N)
int i, j; //счетчики
double* A; //вектор коэффициентов
double Y[5]; //вектор расчетных Y
int main(int argc, char** argv)
{
setlocale(LC_ALL, ".ACP"); //установка поддержки русского языка для вывода на экран
printf("\nВведите порядок регрессионного полинома (1...4) => ");
cin >> m; //вводим порядок регрессионного полинома
A = new double[m + 1]; //выделяем память для вектора А
printf("\nВведите коэффициенты регрессионного полинома:\n");
for (i = 0; i < m + 1; i++) {
printf(" a[%d] = ", i);
cin >> A[i]; //вводим коэффициенты А
}
for (i = 0; i < N; i++) {
Y[i] = A[m]; //Начальное значение переменной Y, задаваемое перед циклом, должно быть равно коэффициенту при Х в старшей степени
for (j = m - 1; j >= 0; j--)
Y[i] = Y[i] * tabl_X[i] + A[j]; //рекуррентно рассчитываем Y
}
printf("\n\nРезультаты расчетов:\n\n");
printf(" X исх. | Y исх. | Y расч.\n");
printf("------------------------------\n");
for (i = 0; i < N; i++)
printf(" %2.0lf | %5.2lf | %5.3lf\n", tabl_X[i], tabl_Y[i], Y[i]);
_getch();
return 0;
}