Алгоритм полных наименьших квадратов в C/C++

Учитывая набор точек P, мне нужно найти линию L, которая наилучшим образом приближает эти точки. Я пытался использовать функцию gsl_fit_linear из научной библиотеки GNU. Однако мой набор данных часто содержит точки, которые имеют линию наилучшего соответствия с неопределенным наклоном (x = c), поэтому gsl_fit_linear возвращает NaN. Насколько я понимаю, для такого рода вещей лучше всего использовать метод наименьших квадратов, потому что он быстрый, надежный и дает уравнение в терминах r и тета (так что x = c все еще можно представить). Кажется, я не могу найти какой-либо код C/C++ для этой проблемы. Кто-нибудь знает библиотеку или что-то, что я могу использовать? Я прочитал несколько исследовательских работ по этому вопросу, но тема все еще немного неясна, поэтому я не чувствую уверенности в реализации своих собственных.

Обновлять:

Я предпринял первую попытку запрограммировать свой собственный броненосец, используя данный код на этой странице в Википедии. Увы, пока безуспешно.

Это то, что у меня есть до сих пор:

void pointsToLine(vector<Point> P)
{
    Row<double> x(P.size());
    Row<double> y(P.size());

    for (int i = 0; i < P.size(); i++)
    {
         x << P[i].x;
         y << P[i].y;
    }

    int m = P.size();
    int n = x.n_cols;

    mat Z = join_rows(x, y);

    mat U;
    vec s;
    mat V;
    svd(U, s, V, Z);

    mat VXY = V(span(0, (n-1)), span(n, (V.n_cols-1)));
    mat VYY = V(span(n, (V.n_rows-1)) , span(n, (V.n_cols-1)));

    mat B = (-1*VXY) / VYY;
    cout << B << endl;
}

вывод из B всегда 0,5504, даже когда мой набор данных изменяется. Также я думал, что на выходе должно быть два значения, поэтому я определенно делаю что-то очень неправильное.

Спасибо!


person McAngus    schedule 22.03.2015    source источник
comment
Это будет закрыто в кратчайшие сроки, так как слишком широкое.   -  person Dirk Eddelbuettel    schedule 22.03.2015
comment
Я не слишком уверен, как еще это выразить. Я ищу реализацию C/C++ общего алгоритма наименьших квадратов для аппроксимации линии. Вот более точное описание того, что это такое.   -  person McAngus    schedule 22.03.2015
comment
См. эту страницу часто задаваемых вопросов по SO и пункт 4: Вопросы с просьбами порекомендовать или найти книгу, инструмент, программную библиотеку. , учебник или другой сторонний ресурс не относятся к теме Stack Overflow, поскольку они, как правило, привлекают самоуверенные ответы и спам. Вместо этого опишите проблему и то, что уже было сделано для ее решения.. В любом случае, на странице Википедии есть код GNU Octave. Теперь используйте, например, Armadillo и напишите его на C++.   -  person Dirk Eddelbuettel    schedule 22.03.2015
comment
Выполните поиск в Интернете по запросу «наименьшие квадраты» или «наименьшие квадраты» без слова «всего». Для линейной аппроксимации подойдет практически любой метод наименьших квадратов.   -  person rcgldr    schedule 23.03.2015
comment
Вопрос поставлен ужасно, но это глупо. ОП не ищет рекомендации, они ищут реализацию или объяснение общих наименьших квадратов в C или C++. Я приземлился здесь в поисках того же самого. Также подход наименьших квадратов не является достаточно хорошим, общий метод наименьших квадратов является допустимым инструментом, когда он подходит к проблеме. Можем ли мы просто быть полезными и правильно ответить на вопрос?   -  person Rehno Lindeque    schedule 30.08.2016


Ответы (1)


Чтобы найти линию, которая минимизирует сумму квадратов (ортогональных) расстояний от линии, вы можете действовать следующим образом:

Линия — это набор точек p+r*t, где p и t — векторы, которые нужно найти, а r меняется вдоль линии. Мы ограничиваем t единичной длиной. Хотя есть и другое, более простое описание в двух измерениях, оно работает с любым измерением.

Шаги

1/ вычислить среднее значение p точек

2/ накапливаем ковариационную матрицу C

    C = Sum{ i | (q[i]-p)*(q[i]-p)' } / N

(где у вас есть N точек и 'обозначает транспонирование)

3/ диагонализовать C и взять в качестве t собственный вектор, соответствующий наибольшему собственному значению.

Все это можно обосновать, исходя из (ортогонального) квадрата расстояния точки q от линии, представленной, как указано выше, которая равна

d2(q) = q'*q - ((q-p)'*t)^2
person dmuir    schedule 23.03.2015