Допустим, перед нами стоит типовая задача регрессии: построить непрерывную зависимость по отдельно заданным точкам. Пускай нам дана зависимость между количеством прочитанных книг по программированию и месячной зарплатой.

Ваня. Прочитал 5 книг, но больше практиковался. Не имеет высшего образования. Сделал свой проект, который стал успешным. Его купила другая крупная фирма. Зарплата руководителя проекта - 300 т.р.

Алексей. Всю жизнь положил на погружение и подготовку. Имеет 2 высших: инженерное и математическое. Прочитал 50 книг. Но устроился работать простым разработчиком встроенных систем в маленькую фирму. Зарплата - 120 т.р.

Вера. С детства увлекается математикой. Закончила мехмат. Прочитала 10 книг. Но открыла в себе талант руководителя. Чаще ставит задачи, иногда тестирует, иногда программирует. Зарплата - 200 т.р.

Ярослав. Гик и людоман. Обожает книги (150) и игры. Когда-то открыл книгу про Ромеро и Кармака. Работал программистом. Сделал инди-игру. Она не выстрелила. Сделал вторую. Она зашла. Работает на основной работе и делает игры дома. Суммарная зарплата - 160 т.р.

Игорь. У него нормально было с математикой в школе. Поддавшись тренду, поступил в универ на ВМСС. Открыл 1 книгу по программированию, но она оказалась неудачной. Решил, что все книги такие. Увлекся КВН, на 2 курсе остался в академе. Бросил универ. Работает в Макдаке. Зарплата - 40 т.р.

Марина. Имеет разряд по шахматам. Прирожденный хакер. Прочитала множество книг (25) по взлому. Устроилась специалистом по цифровой безопасности в крупный банк. Зарплата - 180 т.р.

Феофан. Закончил физтех. Прочитал 30 книг. Устроился специалистом по data science в банк. Зарплата - 80 т.р.

Попытаемся выяснить, как количество прочитанных книг по программированию влияет на зарплату. В задаче будем считать, что месячная зарплата зависит только от этого фактора, хотя в реальности это не так. Упорядочим данные по количеству прочитанных книг.

Изобразим точки нашей дискретной функции на графике в Octave:

x=[1 5 10 25 30 50 150];
y=[40 300 200 180 80 120 160];
plot(x,y,'+');
grid on;
xlabel('x - books count');
ylabel('y - monthly offer');

Линейная регрессия

Мы хотим некое предсказание по "накопленным опытным данным". Какое идеальное количество книг для получения максимальной зарплаты? Построим предсказание в виде линейного графика регрессии. То есть нам надо найти зависимость размера ежемесячной зарплаты y от количества книг в виде . Нам надо найти коэффициенты a и b. В соответствии с выкладками для этого нам надо решить систему линейных уравнений:

Вычислим в Libre Office Calc нужные суммы:
, , ,.

Подставляем полученные значения и получаем систему уравнений относительно a и b:

Выражаем из первого уравнения b, подставляем во второе и находим искомые коэффициенты:

Выразим коэффициенты a и b из системы уравнений в общем виде:

Отсюда получаем формулу для нахождения a и b через суммы:

Напишем программу для нахождения коэффициентов и соберем ее на Visual Studio Community Edition:

#include <vector>
#include <iostream>

#define mypoints std::vector<double>

void computeLinearRegression(mypoints x, mypoints y, double* a, double* b) 
{
    double sum_x = 0, sum_y = 0, sum_xy = 0, sum_x_squared = 0;
    uint32_t num_points = x.size();

    for (uint32_t i = 0; i < num_points; i++)
    {
        sum_x += x[i];
        sum_y += y[i];
        sum_xy += x[i] * y[i];
        sum_x_squared += x[i] * x[i];
    }

    double n = num_points;

    *b = (sum_xy * sum_x - sum_y * sum_x_squared) / (sum_x * sum_x - n * sum_x_squared);

    *a = (sum_y - *b * n) / sum_x;
}

    

int main() {
    mypoints x = {1,  5,   10,  25,  30, 50,  150};
    mypoints y = {40, 300, 200, 180, 80, 120, 160 };

    double slope, intercept;
    computeLinearRegression(x, y, &slope, &intercept);

    std::cout << "Linear Regression Results: ";
    std::cout << "Number of points: " << x.size() << std::endl;
    std::cout << "Slope (a): " << std::fixed << slope << std::endl;
    std::cout << "Y-Intercept (b): " << intercept << std::endl;
    std::cout << "Equation: y = " << slope << "x + " << intercept << std::endl;

    return 0;
}

После запуска получаем вывод:

Linear Regression Results: Number of points: 7
Slope (a): -0.084869
Y-Intercept (b): 157.571343
Equation: y = -0.084869x + 157.571343

Отобразим на графике с точками полученную линейную функцию :

x=[1 5 10 25 30 50 150];
y=[40 300 200 180 80 120 160];
x1=0:0.2:200
y1 = -0.085*x1 + 157.6
plot(x, y, '+', x1, y1);
grid on;
xlabel('x - books count');
ylabel('y - monthly offer');

Мы видим отчетливый тренд на понижение среднемесячной зарплаты, связанный с увеличением прочитанных книг.

Полиномиальная регрессия

Теперь попробуем сделать предсказание в виде степенной зависимости . Нам надо найти коэффициенты a, b и c. В соответствии с методом наименьших квадратов после выкладок имеем следующую систему уравнений:

Вычислим в Libre Office Calc следующие суммы:
, , , ,, ,

Получаем систему уравнений с безумными множителями:

Перепишем систему уравнений в матричном виде и решим ее методом Крамера. Найдем сначала общий определитель, убедимся, что он не равен нулю. Затем вычислим 3 дополнительных определителя, подставляя свободный столбец. Определитель 3 на 3 будем вычислять следующим образом:

На сайте mathprofi.ru есть xls-файл, в котором автоматизирован расчет определителей и нахождения неизвестных методом Крамера. Можно открыть файл в Libre Office Calc и пересохранить в ods и там посчитать все. В Octave определитель матрицы 3 на 3 можно посчитать при помощи функции :

A = [26651 271 7;3543751 26651 271;513711251 3543751 26651]
A =

2.6651e+04 2.7100e+02 7.0000e+00
3.5438e+06 2.6651e+04 2.7100e+02
5.1371e+08 3.5438e+06 2.6651e+04

det(A)
ans = -2.4611e+12

Посчитаем определитель вручную:

Основной определитель не равен нулю. Значит наша система имеет одно решение. Для следующей системы:

Для нахождения неизвестных надо найти общий определитель и дополнительные определители. Например, дополнительный определитель для 1-го столбца выглядит так

Неизвестные найдем по следующим формулам:

Напишем программу для нахождения неизвестных коэффициентов методом Крамера.Так как у нас присутствуют числа порядка (FLT_EPSILON = 2⁻²³ ≈ 1.19e-07 ), надо для вычислений с плавающей запятой использовать тип double (DBL_EPSILON = 2⁻⁵² ≈ 2.20e-16).

#include <vector>
#include <iostream>

#define MATR_SIZE 3
#define SUBMATR_SIZE 2

#define mat_ext std::vector<double>


// квадратная матрица 2 степени, которая получается из исходной квадратной матрицы 3 степени 
// вычеркиванием i-го столбца и j-й строки
mat_ext mysubmatr(SUBMATR_SIZE*SUBMATR_SIZE, 0);

#define SUBMATR(i,j) mysubmatr[i*SUBMATR_SIZE + j]


// временная таблица, в которой i-й столбец заменен столбцом свободных членов
mat_ext tempMatrix(MATR_SIZE*MATR_SIZE, 0);

#define TEMPMATR(i,j) tempMatrix[i*MATR_SIZE + j]

// искомые значения a, b, c в виде вектора
mat_ext solutions(MATR_SIZE, 0);


void makeSubmatrix(mat_ext matr, uint8_t row, uint8_t col)
{
// опредение макроса для удобства работы с линейным вектором в виде матрицы
#define MATR(r,c) matr[(r)*MATR_SIZE + c]

    // сдвиг номера строки подматрицы относительно старой
    int di = 0;

    // сдвиг номера колонки подматрицы относительно старой
    int dj = 0;

    // цикл относительно новой подматрицы
    for (int i = 0; i < SUBMATR_SIZE; i++)
    {
        if (row == i)
            di = 1;

        dj = 0;
        for (int j = 0; j < SUBMATR_SIZE; j++)
        {
            if (j == col)
                dj = 1;

            SUBMATR(i, j) = MATR(i+di, j+dj);
        }
    }
}

// рекурсивное вычисление определителя
double computeDeterminantVector(mat_ext matr, uint8_t n)
{
// опредение макроса для удобства работы с линейным вектором в виде матрицы
#define LOCALMATR(i,j) matr[(i)*n + j]

    if (n == 1)
        return LOCALMATR(0,0);

    if (n == 2)
        return LOCALMATR(0,0) * LOCALMATR(1,1) - LOCALMATR(0,1) * LOCALMATR(1,0);

    double det = 0;
    int sign = 1;

    uint8_t target_row = 0;

    // цикл вычисления детерминанта степени 3 в виде разложения на сумму детерминантов 2 степени
    // раскадываем по первой (нулевой) строке
    for (int j = 0; j < MATR_SIZE; j++)
    {
        // формируем submatrix 
        makeSubmatrix(matr, target_row, j);

        det += sign * LOCALMATR(0,j) * computeDeterminantVector(mysubmatr, n - 1);
        sign = -sign;
    }

    return det;
}

void makeTempMatr(mat_ext mymatr, mat_ext myconst, uint8_t col)
{   
    // опредение макроса для удобства работы с линейным вектором в виде матрицы
#define MYMATR(i,j) mymatr[(i)*MATR_SIZE + j]

    // Copy original matrix
    for (int i = 0; i < MATR_SIZE; i++)
        for (int j = 0; j < MATR_SIZE; j++)
            if (j == col)
                // заменяем i-й столбец свободными членами 
                TEMPMATR(i, col) = myconst[i];
            else
                // другие стобцы просто присваиванием
                TEMPMATR(i, j) = MYMATR(i, j);
}

void solveUsingCramersRuleVector(mat_ext mymatr, mat_ext myconst, uint8_t n)
{
    double mainDeterminant = computeDeterminantVector(mymatr, n);

    // если детерминант равен нулю, тогда > 1 решений 
    if (mainDeterminant == 0)
    {
        std::cout << "The system has no unique solution (determinant is zero).\n";
        return;
    }

    // вычисляем детерминанты
    for (uint8_t i = 0; i < n; i++)
    {
        makeTempMatr(mymatr, myconst, i);

        // неизвестная
        solutions[i] = computeDeterminantVector(tempMatrix, n) / mainDeterminant;
    }

    // вектор найденных неизвестных
    std::cout << "\nSolutions:\n";
    for (int i = 0; i < n; i++) 
    {
        printf("x%d = %.2f\n", i + 1, solutions[i]);
    }
}

int main() 
{
    uint8_t n = MATR_SIZE;

    // 26651.0,     271.0,     7.0,     1080.0,
    // 3543751.0,   26651.0,   271.0,   40440.0,
    // 513711251.0, 3543751.0, 26651.0, 4112040.0

    mat_ext mymatr =  { 26651.0,     271.0,     7.0,
                        3543751.0,   26651.0,   271.0,
                        513711251.0, 3543751.0, 26651.0};

    mat_ext myconst = { 1080.0,
                        40440.0, 
                        4112040.0 };

    solveUsingCramersRuleVector(mymatr, myconst, n);

    return 0;
}

Получаем вывод:
Solutions:
x1 = 0.0098
x2 = -1.6092
x3 = 179.2513

Отобразим график функции

Код для octave:

x=[1 5 10 25 30 50 150];
y=[40 300 200 180 80 120 160];
x1=0:0.2:200;
y1 = 0.0098*x1.^2 - 1.6092*x1 + 179.2513;
plot(x, y, '+', x1, y1);
grid on;
xlabel('x - books count');
ylabel('y - monthly offer');

Полиномиальная регрессия дает совсем другую картину. Мы получаем более интуитивно близкий и согревающий вывод. Чем больше книг прочитано, тем больше зарплата. На этом и остановимся.