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

#include <iostream>
#include <cmath>
#include <iomanip>
#define N 9 //количество експериментов
#define xcp 0 //индекс таблицы в котором храняться средние значения х
#define ycp 1 //индекс таблицы в котором храняться средние значения у
#define tillday 2 //индекс таблицы в котором храняться значения ~y
#define form 3 //индекс таблицы в котором храняться значения выражения
#define a 0 //индекс коэффициента а
#define b 1 //индекс коэффициента b
#define x 0 //индекс коэффициента x
#define y 1 //индекс коэффициента y

using namespace std;

//функция вычисления среднего арифметического елементов массива
double average(double* arr) {
    double res = 0;
    for (int i = 0; i < N; i++)
        res += arr[i];

    return res / N;
}

//функция вычисления среднего геометрического елементов массива
double geometric_mean(double* arr) {
    double res = 1;
    for (int i = 0; i < N; i++)
        res *= arr[i];

    return pow(res, 1./N);
}

//функция вычисления среднего гармонического елементов массива
double harmonic_mean(double* arr) {
    double res = 0;
    for (int i = 0; i < N; i++)
        res += 1/arr[i];

    return N / res;
}

//функция вычисления ~y
double calc_tilday(double* X, double* Y, double _xcp) {
    int tmp_i = 0;
    for (int i = 0; i < N; i++){
        if (_xcp == X[i])
            return Y[i];
        else if (X[i] < _xcp)
            tmp_i = i;
    }

    return ((_xcp - X[tmp_i]) / (X[tmp_i + 1] - X[tmp_i])) * (Y[tmp_i + 1] - Y[tmp_i]) + Y[tmp_i];
}

//функция выведения таблицы на экран
void print_table(double** tabl, int i_func) {
    cout << "Dependency selection table: " << endl;
    cout << "| # |  xcp\t|  ycp\t|  ~y\t|  abs(ycp - ~y) / ~y\t" << endl;
    for (int i = 0; i < N; i++) {
        cout << "| " << i + 1 << ' ';
        for (int j = 0; j < 4; j++)
            cout << "| " << fixed << setprecision(3) << tabl[j][i] << '\t';
        if (i == i_func)
            cout << "+";
        cout << endl;
    }
}

//функция выведения одномерного массива на экран
void print_arr(double* arr, char name_arr) {
    cout << name_arr << " = { " << arr[0];
    for (int i = 1; i < N; i++)
        cout << ", " << arr[i] ;
    cout << " }" << endl;
}

//функция создания таблицы и введения в неё данных
double** build_table(double* X, double* Y) {
    //создаём пустую таблицу
    double** res_tabl = new double* [4];
    for (int i = 0; i < 4; i++)
        res_tabl[i] = new double[N] ();

    //создаём массив врменных значений
    double *tmp_res = new double[3];
    //вычисляем и записываем средние значения Х
    tmp_res[0] = average(X);
    tmp_res[1] = geometric_mean(X);
    tmp_res[2] = harmonic_mean(X);
    //вводим эти значения в таблицу
    for (int i = 0; i < 3; i++)
        for (int j = i; j < N; j += 3)
            res_tabl[xcp][j] = tmp_res[i];

    //вычисляем и записываем средние значения для У
    tmp_res[0] = average(Y);
    tmp_res[1] = geometric_mean(Y);
    tmp_res[2] = harmonic_mean(Y);
    //вводим эти значения в таблицу
    for (int i = 0; i < N; i += 3)
        for (int j = i; j < N; j++)
            res_tabl[ycp][j] = tmp_res[i/3];

    //вычисляем и вводим в таблицу значения ~У и выражения
    for (int i = 0; i < N; i++) {
        res_tabl[tillday][i] = calc_tilday(X, Y, res_tabl[xcp][i]);
        res_tabl[form][i] = fabs((res_tabl[ycp][i] - res_tabl[tillday][i]) / res_tabl[tillday][i]);
    }

    return res_tabl;
}

//функция нахождения индекса минимального елемента
int min(double *arr) {
    int min = 0;
    for (int i = 0; i < N; i++)
        if (arr[min] > arr[i])
            min = i;

    return min;
}

//функция для преобразования елементов массивов Х и У
//для их дальнейшего подставления в формулу нахождения а и b
double** func(double* X, double* Y, int i_func) {
    double** res = new double* [2];
    for (int i = 0; i < 2; i++)
        res[i] = new double[N];
    if(i_func == 0)
        for (int i = 0; i < N; i++) {
            res[y][i] = Y[i];
            res[x][i] = X[i];
        }
    else if (i_func == 1)
        for (int i = 0; i < N; i++) {
            res[y][i] = Y[i];
            res[x][i] = 1. / X[i];
        }
    else if (i_func == 2)
        for (int i = 0; i < N; i++) {
            res[y][i] = Y[i];
            res[x][i] = log(X[i]);
        }
    else if (i_func == 3)
        for (int i = 0; i < N; i++) {
            res[y][i] = log(Y[i]);
            res[x][i] = X[i];
        }
    else if (i_func == 4)
        for (int i = 0; i < N; i++) {
            res[y][i] = log(Y[i]);
            res[x][i] = log(X[i]);
        }
    else if (i_func == 5)
        for (int i = 0; i < N; i++) {
            res[y][i] = log(Y[i]);
            res[x][i] = 1. / X[i];
        }
    else if (i_func == 6)
        for (int i = 0; i < N; i++) {
            res[y][i] = 1. / Y[i];
            res[x][i] = X[i];
        }
    else if (i_func == 7)
        for (int i = 0; i < N; i++) {
            res[y][i] = 1. / Y[i];
            res[x][i] = log(X[i]);
        }
    else if (i_func == 8)
        for (int i = 0; i < N; i++) {
            res[y][i] = 1. / Y[i];
            res[x][i] = 1. / X[i];
        }
    else
        cout << "ERROR: i_func = " << i_func;

    return res;
}

//функция вычисления а и b
double* calc_coef(double* X, double* Y, int i_func) {
    //создаем массив коэффициентов
    double* coef = new double[2];
    //преобразовываем х и у под нужную формулу
    double** newXY = func(X, Y, i_func);

    double sumxy = 0, //сумма х*у
           sumx = 0, //сумма х
           sumy = 0, //сумма у
           sumx2 = 0; //сумма х^2
    //вычисляем суммы
    for (int i = 0; i < N; i++) {
        sumxy += (newXY[x][i] * newXY[y][i]);
        sumx += newXY[x][i];
        sumy += newXY[y][i];
        sumx2 += (newXY[x][i] * newXY[x][i]);
    }

    //вичисляем коэффициенты за формулами(1)
    coef[a] = (sumxy - (1. / N) * sumx * sumy) / (sumx2 - (1. / N) * sumx * sumx);
    coef[b] = (1. / N) * sumy - (coef[a] / N) * sumx;

    //если нужно преобразовываем коэффициенти до нормального вида
    if (i_func == 3) {
        coef[a] = exp(coef[a]);
        coef[b] = exp(coef[b]);
    }
    else if (i_func == 4)
        coef[b] = exp(coef[b]);

    return coef;
}

//функция вывода функции НМК с коэфициентами
void print_func(double* coef, int i_func) {
    cout << "Function MNK: " << endl;
    if (i_func == 0)
        cout << "y = " << coef[b] << " + " << coef[a] << " * x" << endl;
    else if (i_func == 1)
        cout << "y = " << coef[b] << " + " << coef[a] << " / x" << endl;
    else if (i_func == 2)
        cout << "y = " << coef[b] << " + " << coef[a] << " * ln(x)" << endl;
    else if (i_func == 3)
        cout << "y = " << coef[b] << " * " << coef[a] << "^x" << endl;
    else if (i_func == 4)
        cout << "y = " << coef[b] << " * x^" << coef[a] << endl;
    else if (i_func == 5)
        cout << "y = exp(" << coef[b] << " + " << coef[a] << " / x)" << endl;
    else if (i_func == 6)
        cout << "y = 1 / (" << coef[b] << " + " << coef[a] << " * x" << endl;
    else if (i_func == 7)
        cout << "y = 1 / (" << coef[b] << " + " << coef[a] << " * ln(x)" << endl;
    else if (i_func == 8)
        cout << "y = x(" << coef[b] << " + " << coef[a] << " * x)" << endl;
    else
        cout << "ERROR: i_func = " << i_func;
}

int main() {
    //инициализируем значения експерементов
    double* X = new double[N] {23.1, 26.3, 29.8, 31.6, 34.7, 48.1, 55.3, 62.1, 72};
    double* Y = new double[N] {0.99, 0.83, 0.65, 0.59, 0.48, 0.38, 0.27, 0.19, 0.1};

    double **table;
    table = build_table(X, Y);
    print_arr(X, 'X');
    print_arr(Y, 'Y');
    cout << endl;
    int i_func = min(table[3]);
    print_table(table, i_func);

    double* coef = calc_coef(X, Y, i_func);
    cout << endl;
    print_func(coef, i_func);

    return 0;
}