LINUX.ORG.RU

2d фит 2 порядка

 , ,


0

1

Товарищи, ищу что-нибудь уже специфически-оптимизированное под такую задачу:

имеется - координаты множества точек на плоскости, тыщ десять плюс-минус.

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

Gnu scientific library и всякие стандартные фиты в R или очень медленно оптимизируют вблизи оптимальных значений или расходятся в щи при хоть немного косячных стартовых параметров кривой. Я написал кое-как цепочку из разных фитов и отжигов, но это как-то очень всё кустарно выглядит и работает долго. И алхимично.

Было бы круто, если бы оказалось, что у меня корявый гугл и кто-то уже решил такую задачу. Казалось бы, такие задачи должны довольно часто встречаться.

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

ИМХО ты не с той стороны подходишь. Нет такой кнопки «Сделать за<?>бись». Всегда решение будет состоять из кусков.

Теперь оптимизируй своё решение как сумму решений. Больше совета дать не могу, т.к. и так обще рассуждаю. Давай подробности.

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

Ну в случае фитов таких простых функций вполне возможно написать такую кнопку. Для линейных функций, кружков и эллипсов уже написаны. А вот для гипербол и общих 2-порядочных штук почему-то я не нашёл. Будь я нормальным программистом, написал бы либу, но для начала надо хотя бы спросить у нормальных программистов, не велосипед ли это.

И я просто не знаю, какие ещё подробности я могу предоставить. Пример входных данных что ли?

Slesarev
() автор топика
Ответ на: комментарий от ziemin

Не, пространство решений - 5 мерное поле действительных чисел. Ну, скажем, с 4 значащими цифрами будет ок, хотя я считаю float на всякий случай, всё равно сильно не оптимизировать 4 знака.

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

Будь я нормальным программистом

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

Думаю правильным решением будет разбить на подзадачи (не только со стороны программирования, но и со стороны математики).

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

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

Входные данные - массив двумерных целых координат, но их достаточно много, чтобы искать решение на континууме (ну или, технически, на очень плотном счётном множестве)

Slesarev
() автор топика
Ответ на: комментарий от ziemin

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

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

Ты тролишь или тебе правда нужна помощь? Разложи рараболу на одну, основную фигуру и бесконечное множество поправочных.

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

Ну и как мне это поможет?

Вот у меня есть функция нахождения расстояния от точки до кривой, реализованная в виде решения уравнения производной расстояния от точки до точки, лежащей на кривой, и перебора этих решений и в виде перебора расстояний от точки до плотно лежащих точек на кривой. Второй вариант немного тормознее, но для отладки куда удобнее пока.

Метрику я особо не улучшу уже, то есть.

Дальше - фиты в общем виде, они прекрасны и работают отлично в общих случаях и у меня тоже, но часто дают сбои - или промахиваются совсем, или останавливаются в локальных минимумах или просто довольно ровных участках градиента - что офигенно для алгоритмов, справляющихся с общими функциями, но слабовато, когда речь идёт о таких подробно изученных на протяжении, блин, нескольких тысячелетий объектах.

Дальше, если я разобью итоговую кривую на сумму разных объектов, то к этому ансамблю мне придётся один хрен применять фиты в общем виде. Я не понимаю, где профит, объясни, пожалуйста.

Если говорить об изменении вида функции, то можно вообще написать преобразование пространства приводящее кривую в, скажем, прямую, и фитить это преобразование, используя взвешенную метрику расстояния от точек до прямой на образе, но задача останется ровно того же вида.

И да, хуже всего дела обстоят с гиперболами, если есть решение только для них, это уже спасёт дело.

Slesarev
() автор топика

Чет я не понял, тебе надо просто отфитировать множество точек? Ну так это любой матпакет сделать может. Я сам работал только с ROOT, он точно может http://root.cern.ch/root/html/TH1.html#TH1:Fit но я думаю что и остальные не хуже.

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

Дальше, если я разобью итоговую кривую на сумму разных объектов, то к этому ансамблю мне придётся один хрен применять фиты в общем виде.

Я же об это же!

Я не понимаю, где профит, объясни, пожалуйста.

А дело в том, что мелкими можно пренебречь. Если у тебя не дискретная задача и ты её не решаешь в общем виде, то появляется волшебная сигма. Если поправка меньше предела, то её и учитывать не надо!

Главное построить аналитичскую модель поправок к параболе и отбросить малозначащие!

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

Если я всё правильно понял то это функция линейного фита для задач вида min(|f(x)-y|^2), так даже gnuplot умеет

А мне нужно что-то вида min(|f(t)-(x,y)|^2) Или хотя бы min(f(x[7])) при условии, что градиент очень корявый. Может быть, градиент можно исправить другой параметризацией, но я что-то никак не придумаю, как.

Если я неправильно понял, поправьте меня. И в таком случае как называется алгоритм, который там реализован? Может, я уже его пробовал.

Slesarev
() автор топика
Ответ на: комментарий от ziemin

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

Но 4 параметра всё равно подбираются тоскливо долго.

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

данных в вершине гиперболы

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

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

Чувак, сколько раз я должен написать «гипербола» чтобы ты уловил разницу?

Плюс, это только усложнит задачу. Блин, я спрашивал, есть ли готовое решение в либах, трактор, который решит задачу через задницу, я уже написал.

Slesarev
() автор топика
Ответ на: комментарий от ziemin

Ты не строй из себя убогого кадровика мехмата, а поясни, чем замена определяемой 2 параметрами канонической гиперболы определяемым 2 параметрами каноническим эллипсом упростит задачу? (остальные параметры всяко остаются в процедуре афинного преобразования) Даже если оставить в стороне дополнительную операцию по проверке сокращений на малость, я не вижу профита.

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

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

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

Как решают интеграл третьей степени? Через подстановку. I = I^3 + C или как-то так. В общем аналогично решаются все подобные задачи.

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

Это если существует решение для аналогичного объекта более низкого порядка, но кривые второго порядка - это уже ниже некуда просто!

Плюс способ параметризации должен быть найден. У меня была надежда найти его аналитически, но ничего хорошего у меня не вышло. Как я ни пробовал переписать, из всего выходит, что меньше 4 параметров уже не сделать. Более того, приведение в канонический вид - поворот и перенос - 3 параметра, ассимптоты - 1 параметр (угол между ними), реально после поворота и переноса остаётся 1 параметр на фит, но его нужно обрабатывать наравне с поворотом и переносом. Кстати по очереди самосогласовывать я ещё не пробовал, а, возможно, стоило, хотя это ещё большая алхимия, если подумать, некрасиво. Но, в любом случае, ничего проще уже нельзя придумать, эти параметры - реально предельно простые. Единственное, что тут реально есть от гиперболы - это один параметр, меньше некуда. Но вот он-то мне нужен больше всего.

Slesarev
() автор топика
Ответ на: комментарий от morse

Вот и я думаю, что-то тут не так. Кто-то полюбому это делал и не раз. И по-умному.

Это фит гистограммы. Мой объект принципиально отличается тем, что он не однозначная функция - гипербола повёрнута произвольным образом и много где для каждой координаты по одной оси есть две координаты по другой. Плюс, для фита гистограмм обычно используются в качестве метрики только расстояния по одной оси, что проще и быстрее, но я это попробовал и получилось отстойно - хотя я это использовал для автоматического поиска приблизительных стартовых параметров фита, подошло неплохо.

Кстати, привет физикам частиц. Или есть ещё кто-то кто юзает рут? Мне его чуваки из ЦЕРНа показали как Ъ рисовалку графиков, хотя по-моему R и помощнее и коммьюнити у него побольше.

Slesarev
() автор топика

Тред не читал. При фиксированном семействе кривых у тебя линейная задача наименьших квадратов. Скажу по своему опыту: большинство стандартных реализаций(тем более в GSL) плохо будут работать с таким объемом данных: там решается система линейных уравнений, если она у тебя большая то лучше вычислять решения не через стандартное псевдообращение, а просто градиентным спуском или методом сопряженных градиентов. Для начала могу посоветовать зафитить свои данные на подвыборке(где-то точек 200).

Если хочешь, можешь прислать мне данные, попробую посмотреть на них.

И самое главное, тебе нужна именно формула? То есть ты уверен что тебе обязательно нужен параметрическое сглаживание? Может подойдет непараметрический метод? Если да, то опять же, могу показать на твоих данных результат локальной регрессии/сверток/сплайнов.

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

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

И мне обязательно нужно параметрическое сглаживание. Я описываю поведение застрявшей на микрочипе 10нм плёнки жидкой воды метастабильной при любом вакууме и устойчивой к плазме, сжигающей лежащий вокруг вольфрам, мне в приличном журнале за описание неведомой физической хрени сплайном отрежут орган, отвечающий за написание статей, бритвой Оккама, причём тупой стороной.

Вот один набор точек, например:

http://pastebin.com/07zpZuaK

У меня на этом конкретном наборе рассыпался фит почему-то, хотя с такой цепью разных алгоритмов это не удивительно. Зато цепь пофитила на ничотак почти всё остальное и даже чуть меньше половины - на достоверно.

Slesarev
() автор топика

а если так:

[latex]a_1 x^2+a_2 y^2+a_3 xy+a_4 x+a_5 y=1[/latex]

и теперь обычный МНК на него натравить. Ну то есть составить матрицу

[latex]B=(b_{ij}), b_i:=(x_i^2,y_i^2,x_i y_i,x_i,y_i)[/latex]

ну и теперь решить уравнение

[latex]B\underline a=1, \underline a:=(a_1,a_2,a_3,a_4,a_5)^T[/latex]

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

Вот именно так у меня и сделан поиск примерных значений. Работает мимо, увы. Для гипербол с размахом более 90 градусов вообще выдаёт повёрнутую на 90 градусов картинку, что логично, но печально.

Slesarev
() автор топика
Ответ на: комментарий от dikiy

Кстати, почему у меня в браузере эти латехконструкции отображаются в виде их кода а не изящных уравнений?

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

минимизируется по x: |Ax - a|^2_2, но тебе никто не мешает минимизировать по x:

|Ax - a|^2_2 + |Bx - b|^2_2 + |Cx - c|^2_2,

ведь это эквивалентно минимизации

|Fx-f|^2_2,

где F и f --- это вертикально застеченные A,B,C и a,b,c соответсвенно.

Вот один набор точек, например:
http://pastebin.com/07zpZuaK

что ты мне прислал? Там 3 столбца и третий всегда равен 255. Объяни подробнее за свои данные.

И мне обязательно нужно параметрическое сглаживание...мне в приличном журнале за описание неведомой физической хрени сплайном отрежут орган

Я не хочу откланяться от основного вопроса, но разве не все равно, используешь ли ты непараметрический метод, или какую-то непонятную параметрическую функцию, не имеющую отношения к физике задачи?

maggotroot
()
Ответ на: комментарий от Slesarev

Кстати, почему у меня в браузере эти латехконструкции отображаются в виде их кода а не изящных уравнений?

установи плагин от Эдика.

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

А не забыл ли ты константый член? Там получится 6 неизвестных. Прошу заметить, что тут ты не можешь просто так поделить на коэф. при единичке.

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

Чтобы так сделать, мне придётся разделить переменные, а в случае с кривыми 2 порядка это чревато проблемами, разве нет? Что-то я никак не соображу, как это сделать грамотно.

Третий столбец это шлак, выхлоп предыдущего шага, сделанного сторонним софтом. Мне проще его игнорить, чем что-то делать с предыдущим шагом. Прошу прощения, я к нему так привык уже, что не заметил.

А про параметрическую функцию - я привязываю все 2 параметра этой гиперболы к настоящей физике и объясняю, почему это именно гипербола или не гипербола, так что они мне нужны. Ассимптоты реально существуют как геометрические объекты, для отклонения вершины гиперболы от их пересечения у меня есть потенциал, который я пытаюсь вычислить и соотнести поведение результата с теорией.

Slesarev
() автор топика
Ответ на: комментарий от maggotroot

Если кривая не вырождена, то могу. А она не вырождена, иначе Ландау в гробу перевернётся.

Slesarev
() автор топика
Ответ на: комментарий от dikiy

А я смогу после этого мануалы по латеху подглядывать в интернетах или все теги отрендерятся сразу?

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

Пожалуйста, пришли мне данные. Потому что я не пойму что от чего тебе надо строить. Где у тебя 2d фит(должно быть 2 признака и 1 регрессор на точку)? Где у тебя фит к разным «осям»? Что значит что надо разделить переменные, это типа ты говоришь, что часть переменных описывается кривой с одними коэффициентами, а часть точек уривой с другими коэф?

Объясни нам кратко и содержательно, какие у тебя данные, и что ты хочешь получить. И желательно выложи данные.

maggotroot
()
Ответ на: комментарий от Slesarev

А я смогу после этого мануалы по латеху подглядывать в интернетах или все теги отрендерятся сразу?

я не распарсил.

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

a=
   -3.4044e-06
    2.3730e-06
   -2.5611e-07
    4.4723e-03
   -2.0253e-03

с формулой: [latex] a_1 x^2 + a_2 y^2 + a_3 xy + a_4 x + a_5 y = 1[/latex]

со среднеквадратичным отклонением 1.88.

dikiy ★★☆☆☆
()
Последнее исправление: dikiy (всего исправлений: 1)
Ответ на: комментарий от maggotroot

Про 2d фит я имел в виду то, что у меня не 2мерное пространство признаков и 1 регрессор, а то, что у меня 2мерное пространство данных и надо найти кривую, которая их аппроксимирует. Функции в явном виде я не имею, иначе бы я точно поюзал бы МНК-рутины.

Формат данных:

в каждой строке:

x y z

(x,y) - координаты точки

z - левый шлак, игнорировать

нужно найти такие a, b, что после поворота и сдвига картинки точки лежали наиболее близко к гиперболе, например, такого вида:

x=+-a*ch(t)

y=+-b*sh(t)

t=[-Inf,+Inf]

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

О, вот это круто, а я её не видал ни разу. Уж не знаю, решит ли это мою задачу, но пару годных оптимизаций я там уже вижу. Дочитаю - отпишусь.

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

В этой пдфке рассмотренно то что тебе предложил dikiy, просто вместо квадрата l_2 нормы минимизируется l_2 норма, тем самым ты получаешь нелинейную оптимизацию.

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

Нет, там ещё есть градиентное взвешивание. Я эти веса пытался сам найти, но это долго и трудно - определить их и показать хотя бы себе что они корректные. А квадрат или не квадрат нормы минимизировать это фигня, и поменять не сложно. Ещё вопрос что корректнее, кстати.

UPD О, я не внимателен, там таки минимизация квадратов рассмотрена.

И ещё там другие методы, я их тоже попробую, посмотрю, какой самый приятный.

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

Это градиентное взвешивание не что иное как IRLS. Прсто метод для минимизации такой нелинейной байды, где на каждом шаге мнк решается. Сейчас тебе скину что наваял на матлабе.

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

Да, я вижу, что это итеративная шняга. Там на следующей странице неитеративный грубоватый вариант дан. Годные решения вне аппендикса. Но всё равно присылай, только плейнтекстом код, самого матлаба у меня нету.

Господи,

The first thorough investigation of the ellipse fitting problem was done in the middle 1990’s

Чем эти математики занимаются, если этой проблеме меньше 20 лет?

Slesarev
() автор топика
Ответ на: комментарий от Slesarev

Блииин, я затупил, надо было сразу данные построить, сразу все стало бы понятно. Короче у тебя подгонка неявно заданной кривой, ясен-красен, что это нелинейная оптимизация! Ты можешь это решать много чем, но я использовал метод левенберга-маркудта, который и создан для нелинейных наименьших квадратов.

В общем тут тебе график с моими коэф и невязкой https://www.dropbox.com/sh/ogc6pyppq5p7co9/Xva-WU6VlQ (пнгшка и епс с одним и тем же). Правда я немного твои данные подстриг: там куча подряд идущих одинаковых значений второй колонки, которым соотв. разные значения первой колонки, я оставил по одному значению, усреднив второе. Ну думаю идею ты понял. Код я сейчас не смогу выложить: у меня там свои матлабовские утилитки которые помогают всю эту нелинейную фигню решать, а я сейчас в отпуске. Если хочешь, я могу в теченеие трех дней тебе на питоне с сайпи переписать. Если тебе срочно, то могу постараться к завтра сделать.

Дай нать, то ли это чты ты хотел.

Пс про ту пдфку: там есть кое-что полезное, но этот азиат(как обычно азиаты и делают) наводит ненужное наукообразие в элементарных вещах.

The first thorough investigation of the ellipse fitting problem was done in the middle 1990’s

Ну он видимо идиот, этим еще Кеплер занимался.

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

Извини что запутываю, это не нелинейная задача. Ну точнее она может стать нелинейной, если плохо ввести ограничения. Короче читай тут, здесь хорошо описано http://research.microsoft.com/en-us/um/people/zhang/Papers/ZhangIVC-97-01.pdf

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

Азиат там сеньёр, а первым автором некий господин Чернов. Конечно, статья написана довольно тоскливым языком, но она походу не в журнал какой-то а для себя, и если считать её литобзором с выписанными самыми нужными формулами, то оно очень ок. Например, алгоритм поиска расстояния до кривой описанный неким Эберли, люто обгоняет всё, что я пробовал, судя по всему, и это приятно. И формулы 6-8 это что-то великолепное просто, чуваки наверняка потратили изрядно времени и травы, чтобы так всё элегантно выписать.

А я про Кепплера тоже подумал, но он всё же наверное не доводил дело до алгоритмов. Небось даже Коперник что-то по теме сделал.

Я вот пишу тихонечко реализацию их ерунды (мне не к спеху, вряд ли кто-то ещё достаточно упорот, чтобы обнаружить и опубликовать эффект, который я обсчитываю - так что я тебе отпишусь когда из отпуска вернёшься, если что-то у меня не срастётся) и думаю - как же весь 20 век астрономы без таких алгоритмов жили? Нет, что-то тут не то. Хотя астрономам-то небось пофигу, если параметры будут определены с ошибкой и это займёт годы процессорного времени, то потом их можно будет тихонько подгонять и подгонять тупо закидывая новыми данными, да и орбиты там довольно сильно от этих Кепплеровских фигур отличаются.

Попутно вопрос - модифицированный Ньютон-Гаусс упоминается в статье и я к нему привык, а вот сейчас погуглил и обнаружил, что он устарел уже по мнению некоторых полвека назад и щас рулят всякие какие-то Nelder-Mead и Broyden-Fletcher-Goldfarb-Shanno и simulated annealing (хотя первый, я так понял, это для детей и градиенты получает численно). Левенберг-маркудт правда уже не лучший, или это фигня всё? (надо бы почитать, пожалуй, что это за алгоритмы, наконец)

Slesarev
() автор топика
Ответ на: комментарий от maggotroot

Коллега, а оставьте мне свою почту, я к вам с дурацкими вопросами буду приставать, если что. Моя - как моё имя тут только последняя буква дубльвэ и на гмейле.

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