Пример метода наименьших квадратов c использованием библиотеки gsl

Question

Прошу поделиться примером МНК для полинома n степени с использованием библиотеки gsl. Я смотрел в документации "least square method", но поиск не дал результата. Видать как-то по другому называется данный метод в документе.

Answer 1

Я посмотрел свои исходники, но в них я вручную составляю систему уравнений для МНК, а библиотекой Eigen её решаю. Возможно, в ней есть и более подходящие инструменты без лишних строк кода, но мой код такой:

double x0 = model->rowCount();
double x1 = 0, x2 = 0, x3 = 0, x4 = 0, x5 = 0, x6 = 0;
double x0y = 0, x1y = 0, x2y = 0, x3y = 0;
for (int i = 0; i < x0; i++)
{
    x1 += qPow(model->calibPoint(i).Frequency, 1);
    x2 += qPow(model->calibPoint(i).Frequency, 2);
    x3 += qPow(model->calibPoint(i).Frequency, 3);
    x4 += qPow(model->calibPoint(i).Frequency, 4);
    x5 += qPow(model->calibPoint(i).Frequency, 5);
    x6 += qPow(model->calibPoint(i).Frequency, 6);

    x0y += model->calibPoint(i).Pressure;
    x1y += qPow(model->calibPoint(i).Frequency, 1) * model->calibPoint(i).Pressure;
    x2y += qPow(model->calibPoint(i).Frequency, 2) * model->calibPoint(i).Pressure;
    x3y += qPow(model->calibPoint(i).Frequency, 3) * model->calibPoint(i).Pressure;
}
Eigen::Matrix<double, 3, 3> eqLeft;
Eigen::Matrix<double, 3, 1> eqRight;

eqLeft << x0, x1, x2,
          x1, x2, x3,
          x2, x3, x4;

eqRight << x0y, x1y, x2y;

Eigen::Matrix<double, 3, 1> solution = eqLeft.partialPivLu().solve(eqRight);

//Существует ли решение системы уравнений
double relative_error = (eqLeft * solution - eqRight).norm() / eqRight.norm();
if (relative_error > 1e-9)
{
    //emit SomeError;
    return;
}
if (std::isnan(relative_error))
{
    //emit SomeOtherError;
    return;
}

//Теперь решение можно получить через solution(0), solution(1), solution(2)

БЛОГ НА HUSL

Пример метода наименьших квадратов c использованием библиотеки gsl

Ответы (1 шт):