LINUX.ORG.RU

Расчет ошибки в 2х параметрическом линейном МНК без матриц.

 ,


0

1

Вот есть расписанный 2х параметрический линейный МНК без матриц: https://math.stackexchange.com/questions/1657030/fit-plane-to-3d-data-using-least-squares

Но там нет выражения для расчета квадратичной ошибки. Я попытался это сам подумать, но это вышем моих способностей и доступного времени. И там добавить вопрос в комментарий не могу - не хватает репутаций.

Может кто знает, как это можно посчитать?

Я подобный вопрос про этот мнк здесь уже спрашивал как-то, в результате нашлась та ссылка. Но расчет ошибки сам так и не осилил.

Но там нет выражения для расчета квадратичной ошибки.

Для параметров модели? Я предположу, что для параметров модели.

Может кто знает, как это можно посчитать?

Без матриц? Не советую.

Пусть координаты точек лежат в массиве X[i,j] (i-я точка, j принадлежит {1, 2}), а линейная модель должна предсказать третью координату y[i], оценив параметры a[j]: y[i] = a[1] * X[i,1] + a[1] * X[i,2] (все координаты центрированы, чтобы не заботиться о постоянном слагаемом).

Вам нужно взять оценку стандартной ошибки третьей координаты (либо ошибку измерения, если она есть, либо модуль остатка, если её нет) σ(y[i]), взять выражения для каждого параметра a[j] и посчитать выражения с частными производными a[j] по y[i]:

σ²(a[j]) = ∑[i] σ(y[i])² (∂a[j]/∂y[i])²

(Это формула 15.4.12 из третьего издания Numerical Recipes.)

С матрицами получается очень легко, достаточно взять C = (X'X)^-1 и получить σ²(a[j]) = C[j,j] (15.4.15), а без матриц придётся расписывать частные производные из уже многоэтажных формул по Вашей ссылке. Если лень делать это вручную, возьмите Maxima (с дифференцированием она должна справиться). На C++ с перегрузкой операторов можно также сделать дуальные числа, но это для эстетов.

AITap ★★★★★
()
Ответ на: комментарий от AITap

Ту ошибку, что «residual sum of squares (RSS)»

То, что с матрицами это легче, это понятно. А без матриц следующие плюсы:

1. Производительность и расчет в один проход без необходимости сохранения матриц.

2. Слагаемость результата - все компоненты представляют собой суммы, и можно взять два набора результатов расчета от двух разных данных и сложить между собой. И вычитаемость. Если есть расчет больших данных, и из пачки делается исключение небольшого набора, то пересчитывать большие данные не обязательно, а можно вычесть расчет исключаемой пачки.

Последний пункт очень удобен для расчетов классификаций - деление множеств данных на пачки.

Для однопараметрического уравнения у меня такие расчеты:

// добавление значения в расчет, sg это знак добавления,
//   и если -1 значит вычитание.
void Corr::add_value(double x, double y, double w, double sg) {
    su_w += w * sg;
    su_w_sqr += w * w * sg;
    su_x_y += x * y * w * sg;
    su_x += x * w * sg;
    su_y += y * w * sg;
    su_x_sqr += x * x * w * sg;
    su_y_sqr += y * y * w * sg;
}

// сложение и вычитание расчитанных пачек 
void Corr::add_other(const Corr &o, double k_other) {
    su_w += o.su_w * k_other;
    su_w_sqr += o.su_w_sqr * k_other * k_other;
    su_x_y += o.su_x_y * k_other;
    su_x += o.su_x * k_other;
    su_y += o.su_y * k_other;
    su_x_sqr += o.su_x_sqr * k_other;
    su_y_sqr += o.su_y_sqr * k_other;
}

// это расчет коэффициентов
PairDbl Corr::get_a_b() const {

    assert( x_more_zero(su_w) );

    double divider = su_w * su_x_sqr - su_x * su_x;

    if (is_zero(divider))
        return {0, su_y / su_w};

    double d1 = su_w * su_x_y - su_x * su_y;

    double a = d1 / divider;
    double b = (su_y - a * su_x) / su_w;

    return {a,b};
}

// а это расчет ошибки.

    auto [a,b] = get_a_b();

    double d3 =
         su_y_sqr
         - 2. * a * su_x_y
         - 2. * b * su_y
         + 2. * a * b * su_x
         + a * a * su_x_sqr
         + b * b * su_w;

    double variance = d3 / su_w;
    
    double error = std::sqrt(variance);

Хочу в подобном стиле сделать для двухпараметрического. Без расчета ошибки оно не имеет смысла начинать.

victor79
() автор топика
Последнее исправление: victor79 (всего исправлений: 2)
Ответ на: комментарий от victor79

Ту ошибку, что «residual sum of squares (RSS)»

А, сумму квадратов остатков, ∑(ŷ-y)². Чтобы её вывести, достаточно подставить выражение для ŷ и раскрыть скобки в квадрате суммы:

∑[i] (yhat[i] - y[i])^2 =
= ∑[i] (X[i,1] * a[1] + X[i,2] * a[2] - y[i])^2 =
= ∑[i] (X[i,1] * a[1] + X[i,2] * a[2] - y[i]) * (X[i,1] * a[1] + X[i,2] * a[2] - y[i]) =
= ∑[i] {
   X[i,1] * a[1] * X[i,1] * a[1] + X[i,1] * a[1] * X[i,2] * a[2] - X[i,1] * a[1] * y[i] +
 + X[i,2] * a[2] * X[i,1] * a[1] + X[i,2] * a[2] * X[i,2] * a[2] - X[i,2] * a[2] * y[i] -
 - y[i] * X[i,1] * a[1] - y[i] * X[i,2] * a[2] + y[i] * y[i]
} = ∑[i] {
 X[i,1]^2 * a[1]^2 + X[i,2]^2 * a[2]^2 + y[i]^2 + 2 (
  X[i,1] * a[1] * X[i,2] * a[2] - y[i] * X[i,1] * a[1] - X[i,2] * a[2] * y[i]
 )
}

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

AITap ★★★★★
()
Ответ на: комментарий от AITap

Супер! Спасибо! Понял, и понял принцип получения. Думал производные нужно выводить, которые я уже совсем забыл. А сами коэффициенты в той ссылке расписаны. Теперь сделаю наконец таки этот класс, давно пытался подступиться.

victor79
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.