Среднеквадратичное приближение функций

Секретный ресурс:

http://www.gnu.org/software/gsl/manual/html_node/Interpolation.html

так это вроде не интерполяция ? по идее вот это больше походит ? http://www.gnu.org/software/gsl/manual/html_node/Least_002dSquares-Fitting.html

там нечего изобретать. Все что тебе надо - решенить систему уравнений.

допустим будешь приближать полиномами n степени:

дано k точек в виде вектора x и вектора y, и пусть k>n.

берем матрицу A_{ij}=(x_i)^j, i=1..k, j=0..n

имеем систему Ac=y, где c - вектор коэффициентов полинома.

Очевидно, что она в общем случае не имеет решений. Но имеет наилучшее квадратическое приближение, которое можно посчитать решив уравнение

A^t*A*c=A^t*y, где A^t - транспонированная матрица.

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

Что-то похожее было в курсе Machine Learning. Там обычная линейная регрессия.

я думаю, что в gsl можно инвертировать матрицы и умножать их? :)

оно ? http://www.gnu.org/software/gsl/manual/html_node/QR-Decomposition-with-Column...

http://reference.wolfram.com/mathematica/ref/Fit.html

Да элементарно все: сначала в Octave подбираешь дельный алгоритм. Потом его переносишь в сишный код. В GSL есть функции для аппроксимации методом наименьших квадратов. Я в своей fitsview так гауссианки вписывал. Хотел так же сделать для окружностей, но там проблема — две точки с бесконечной производной.

оно ? http://www.gnu.org/software/gsl/manual/html_node/QR-Decomposition-with-Column...

не совсем. Ищи просто умножение матриц. Когда ты получишь матрицу A: A_{ij}=(x_i)^j, то вектор коэффициентов будет

c=inverse(A^t*A)*A^t*y

Используй rloess в матлабе или R. Реализация rloess на си с gsl у меня есть, если надо. Короче локальную регрессию.

Вспомнил (точнее, подсмотрел): gsl_multifit — итеративный решатель функций с несколькими параметрами.

Изобретать велосипед не хочется ( как сделал автор статьи )

«Велосипед» — это один цикл и несколько арифметических операций? OMFG...

МНК на программируемых калькуляторах и то требует несколько десятков байт только...

у меня лабораторка такая была в университете, все элементарно там

МНК на программируемых калькуляторах и то требует несколько десятков байт только...

Разве что для линейной аппроксимации (или аппроксимации зависимости, которую можно свести к линейной).

Оно элементарно все конечно, но не хочется заново все переписывать если в будущем будет изменем метод/алгоритм решения. Так что решение действительно оказалось элементарно: в gsl в примерах есть файл fitting2.c, где как раз эта задача и приводится. Достаточно ввести массив значений на вход программы - на выходе получаем искомую кривую. Там используется gsl_multifit_wlinear.

matlab: polyfit(x,y,n)

gsl_multifit_wlinear

Сдается мне, что это — кусочно-линейная аппроксимация. А в ТЗ была аппроксимация функцией.

Разве что для линейной аппроксимации (или аппроксимации зависимости, которую можно свести к линейной).

Угу. Для любой представимой двухпараметрической. Через преобразование f(y) = a+b*f(y). Путём подбора разных простых f(x) (1/x, x², sqrt(x), log(x) и т.д.) можно находить готовые простые функции, весьма точно описывающие большинство осмысленных кривых. На МК-61 у меня была программа, таким образом автоматом перебирающая 16 функций (4 функции по x и y) и я любил развлекаться, подбирая зависимости ко всяким табличным данным, типа зависимости давления насыщенного водяного пара от температуры и т.п. Точность была весьма большая.

Ну а если тупо многочленами, то ресурсов калькулятора хватало, наверное, до 5..7 порядка, где-то. Но я никогда не любил такую аппроксимацию, больно у неё фальсифицируемость низкая :)

Там и посложнее можно. Я помню, на лабораторках по физике мы делали аппроксимацию всяких функций вида y = a*exp(bx), y = a*ln(bx) и т.п. — путем несложных преобразований они легко линеаризуются, затем находятся коэффициенты, преобразуются в нужный вид и получается зависимость.

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

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

y = a*exp(bx)

Угу, это же log(y) = a' + bx :)

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

Но в общем случае придется решать оптимизационную задачу :)

Берём матрицу преобразований по X и Y и ранжируем результат по коэффициенту корреляции или по дисперсии (зависит от исходных данных).

Разве что для линейной аппроксимации (или аппроксимации зависимости, которую можно свести к линейной).

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

Ну или сплайнами сгладить.

сплайнами ты интерполяцию сделаешь. Это совсем другое.

только что глянул в GSL. ТСу gsl_multifit_linear доктор прописал.

Не, можно и аппроксимацию сплайнами сделать.

С gsl_multifit_fdfsolver интереснее же!

Не, можно и аппроксимацию сплайнами сделать.

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

С gsl_multifit_fdfsolver интереснее же!

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

То есть лучше все таки стараться линейными комбинациями приближать. Или хотя бы стартовую позицию с умом выбирать.

да и вообще, имхо gsl не нужен в 99% случаев :) Его с успехом заменяет octave.

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

Использовать такое надо с умом

Это точно: я, наверное, дня полтора убил, пока научился гауссианку вписывать. Причем, первые часа три тупо втыкал, почему результат расходится, пока не заметил, что производную неправильно записал. Зато получаем хорошую (если правильно выбрать начальные условия, как ты справедливо заметил) аппроксимацию функции >2 параметров.

gsl не нужен в 99% случаев :) Его с успехом заменяет octave

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

гугл, второй результат сверху - то что тебе надо.

cubic splines orthonormal

а блин. Ты про готовое решение говорил.

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

И очень странно, что до сих пор никто этого не сделал: задача-то очень важная. Астрофизики уже лет 20 (или сколько там уже ПЗС-матрицам) страдают, а все равно пока единственный надежный способ — вручную натыкать опорных точек и проводить кривую по ним. А хочется-то автоматизировать все это.

// я, кстати, один метод снижения шумов пытался как-то реализовать, но пока ресурсов недостаточно: надо хотя бы ~64ГБ оперативки, иначе считается две недели.

хм. которкое гугленье вывело меня на NAG fortran library. Нужная тебе функция называется e02baf - аппроксимация кубическими сплайнами.

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

Моя жопа говорит мне, что данная проблема является выпуклой -> можно ЧМО использовать для оптимизации параметров.

а уж медианная так вообще темный лес!

что такое медианная?

можно ЧМО использовать для оптимизации параметров

Какое ЧМО?

Статьи, кстати, на эту тему есть. Но там жуткие уравнения. Их пока в нормальный вид приведешь, рехнуться можно!

что такое медианная?

Вычисляем медиану отклонения экспериментальных данных от аппроксимирующей функции. Оптимальными параметрами считаем те, у которых эта медиана будет минимальной.

можно ЧМО использовать для оптимизации параметров

Какое ЧМО?

численные методы оптимизации. Я немного неточно выразился :) Имел в виду к примеру банальный градиентный спуск.

Статьи, кстати, на эту тему есть. Но там жуткие уравнения. Их пока в нормальный вид приведешь, рехнуться можно!

математики вам там нужны? :)

Вычисляем медиану отклонения экспериментальных данных от аппроксимирующей функции. Оптимальными параметрами считаем те, у которых эта медиана будет минимальной.

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

А, понял. Ну, как появится нужда — почитаю еще, как оптимизировать задачу.

математики вам там нужны? :)

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

походу такая задача вообще не будет иметь единственного решения

Естественно.

А то, здесь физиков — пруд пруди, а вот математиков хороших хрен найдешь.

хых. как закончу учиться - пошлю свое резюме :)

походу такая задача вообще не будет иметь единственного решения

Естественно.

а какой тогда от нее толк?

как закончу учиться - пошлю свое резюме

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

а какой тогда от нее толк?

От решения некорректных задач тоже толк бывает.

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

а вот кстати интересно. Мне как раз надо будет практику пройти через пол-года где-то.

Вообще в принципе это реально к вам на практику приехать? Ну и если реально - сколько примерно стоит добраться из Европы?

Вообще в принципе это реально к вам на практику приехать?

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

сколько примерно стоит добраться из Европы?

Сколько до Москвы или МинВод — не знаю. А от Москвы самолет стоит тысяч 5, поезд — от 1.5т.р. (смотря какой поезд и какой класс). Из МинВод на автобусах (весь день займет и рублей 400), на такси (около 2.5т.р.) или, если получится, кто-нибудь подберет (если будет оказия с дежуркой, то вообще бесплатно, а так — от 1.5 до тех же 2.5т.р.).