Подгонка функции (curve fitting)

[Евгений Машеров]

Матрица одна. Эн на Эн размерностью. На каждом шаге обращается она одна.

[ewert]

Mavlik в сообщении #1044461 писал(а):

$\alpha=\frac{1}{_{t_{1}}}\cdot \ln(\frac{sz}{\sin (\omega\cdot t_{1})})$

Это не есть хорошо -- возможны проблемы с точностью при переходах через ноль. И вылет запросто возможен из-за случайной отрицательности в окрестностях нулей (частота/фаза-то определены неточно).

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

++++++++++-+------------...

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

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

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

Насчёт программы. Судя по типу single Вы явно используете Borland/Turbo Pascal и потому испытываете нехватку памяти. Выделяйте её в куче или переходите на Free Pascal; он хоть и чуть менее удобен, но зато проблемы с памятью снимает полностью.

[Deggial]

Евгений Машеров в сообщении #1044570 писал(а):

Матрица одна. nxn.

Евгений Машеров, замечание за хроническое неоформление формул. В следующий раз буду сносить ваши посты в Карантин.

[Евгений Машеров]

А Вы уверены, что это формула? $N=n\cdot n$ это формула. А здесь лишь указание на размерность.

-- 12 авг 2015, 22:47 --

Тип single, 32-битный плавающий, вообще стоит избегать. Затрудняюсь сказать, где он может быть полезен - для хранения больших массивов при острой нехватке памяти и отсутствии больших вычислений с ними? На старых машинах использовалась разрядность 36 бит, и это был minimum minimorum, а то и 48 бит. Но с полвека назад скрестилиежа с ужоммашины для экономических расчётов с байтовой организацией и числогрызы с плавающей точкой, вот и появились недостаточный по точности single и избыточный double/

[ewert]

Евгений Машеров в сообщении #1044858 писал(а):

А здесь лишь указание на размерность.

И читается это указание всеми ежами исключительно как "нхн". А если б Вы были не автором, а ежом -- Вам бы самому такое прочтение понравилось?...

Евгений Машеров в сообщении #1044858 писал(а):

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

Я же сказал, где -- в Borland/Turbo Pascal. ТС захотелось пять тыщ отсчётов, а они умещаются в разрешённые 64 килобайта только при типе single. Ни при real, ни при double, ни тем паче при extended они в стандартную 286-ю память не умещаются.

[Евгений Машеров]

(Оффтоп)

Если бы я был ежом - я бы начал с выяснения вопроса, не вкусный ли это червячок? Нет? Ну хоть что-то съедобное? Тоже нет? Тема бесперспективная и абсолютно недиссертабельная.
Увы, концептуальная парадигма ежей несколько прагматична...

Давненько я не брал в руки Turbo Pascal. Хотя не так давно, как могло быть. В обозримом прошлом пришлось корректировать софт к одному старому прибору, где разработчик кое-что наврал. Любезное предоставленные им исходники компилировались только на TP, более того, исключительно на TP 5.5 (что-то он там на ассемблере вставлял хитрое). Но даже на TP, помнится, были фокусы, позволявшие брать из "кучи" более 64Кб одним куском, и даже брать из extended memory, до 4 Мб (не уверен, но возможно, его ASM-вставки как раз для этого и были, хотя, может, для фокусов с видеопамятью). Впрочем, более осмысленно идти по стопам Древних Титанов, храня такие, не влезающие в ОП, данные, во внешней памяти и подгружая по надобности (Великие Древние ухитрялись обращать матрицу размерности тысяча на тысяча при 4096 словах памяти...). И даже если хранить исходные данные в single (а они, подозреваю, снимаются с АЦП, где разрядность 12 или 14, может, 16, и можно хранить, как целочисленные отсчёты по 2 байта), то вычислять лучше всё же в double. Обращая матрицу, в частности. Маленькие ошибки могут там изрядно вывернуть обратную матрицу.

Возвращаясь к Левенбергу-Марквардту. Прежде всего - это приближение к методу Ньютона. Но упрощённое. Не требующее вычисления вторых производных от приближающей функции. Матрица вторых производных от квадратичной целевой функции (суммы квадратов невязок) включает первые производные от приближающей функции и вторые производные от неё. Упрощение в Л-М состоит в том, что вторые производные мы волевым решением приравниваем к нулю, избавляясь от $n^2$ вычислений вторых производных на каждой итерации и такого же количества выписываний аналитических выражений для них. Обращаемая матрица при этом может быть вырождена или близка к вырожденной, и мы ещё более загрубляем, прибавляя к её диагонали положительные числа (одинаковые или пропорциональные диагональным элементам).
То есть нам надо уметь считать $n$ производных по параметрам.
Если функция имеет вид $y=e^{-at}(I_0\sin\omega t+I_1\cos\omega t)$
(я настаиваю на необходимости учёта фазового сдвига - даже если Вы уверены, что при $t=0$ непременно $y=0$, ничто не гарантирует, что Ваш первый отсчёт, снятый с АЦП, будет соответствовать $t=0$; кроме того, использование представления в виде суммы синуса и косинуса, синфазной и квадратурной составляющей, несколько упрощает работу сравнительно с $\sin(\omega t+\varphi)$)
то производные будут
$\frac{\partial y}{\partial a}=-t e^{-at}(I_0\sin\omega t+I_1\cos\omega t)$
$\frac{\partial y}{\partial \omega}=e^{-at}(I_0 t \cos\omega t-I_1 t \sin\omega t)$
$\frac{\partial y}{\partial I_0}= e^{-at}\sin\omega t$
$\frac{\partial y}{\partial I_1}= e^{-at}\cos\omega t$
Имея какие-то начальные приближения для параметров (лично я бы в Вашем случае частоту оценил бы через число максимумов, минимумов или пересечений ноля, а декремент затухания - проведя линейную регрессию через логарифмы максимумов и/или абсолютных величин минимумов, наклон даст Вам параметр $a$, а свободный член будет логарифмом $I$), вычисляем значения функции в точках при этих параметрах и их отклонения от наблюдаемых, составляющие вектор $u$, а значения производных в этих точках - матрицу Z (сержант-вычислитель Равненько из батареи к-на Очевидность подсказывает, что в программе их хранить необязательно, особенно если у Вас дефицит памяти, о векторе и матрице речь идёт лишь для удобства изложения, на практике лучше держать только матрицу $R=Z^TZ$ и вектор $w=Z^Tu$, существенно меньшей размерности, накапливая в них значения, двигаясь по $t$)
Тогда поправка к значениям параметров получается, как $\Delta=(Z^TZ)^{-1}Z^Tu=R^{-1}w$, но тут может быть вычислительная проблема, если матрица $R$ будет близка к вырожденной, поэтому пользуются более устойчивой вычислительно поправкой $\Delta=(Z^TZ+D)^{-1}Z^Tu=(R+D)^{-1}w$, где $D$ - диагональная матрица, элементы которой либо равные малые числа, либо умноженные на малое число диагональные элементы матрицы $R$
(насколько малые - ну, вот я увеличивал диагональные на один процент и сверх того прибавлял одну миллионную).
Скорректировав значения параметров - повторяем расчёт, следя за сходимостью.

[worm2]

Евгений Машеров в сообщении #1044858 писал(а):

Тип single, 32-битный плавающий, вообще стоит избегать. Затрудняюсь сказать, где он может быть полезен

Побойтесь ктулху! Он отлично подходит почти для всех исходных данных. Редко где измерительные приборы могут обеспечить точность больше 7 значащих цифр. Да и простые операции с исходными данными (2-3 штуки) чаще всего не приводят к неприемлемому уменьшению точности. Я с геофизическими данными работаю, где почти всё в single. И даже там, где нужно получить разность двух близких измерений (частая операция), single всё ещё отлично работает.

[Евгений Машеров]

worm2 в сообщении #1044946 писал(а):

Евгений Машеров в сообщении #1044858 писал(а):

Тип single, 32-битный плавающий, вообще стоит избегать. Затрудняюсь сказать, где он может быть полезен

Побойтесь ктулху! Он отлично подходит почти для всех исходных данных. Редко где измерительные приборы могут обеспечить точность больше 7 значащих цифр. Да и простые операции с исходными данными (2-3 штуки) чаще всего не приводят к неприемлемому уменьшению точности. Я с геофизическими данными работаю, где почти всё в single. И даже там, где нужно получить разность двух близких измерений (частая операция), single всё ещё отлично работает.

(Оффтоп)

А чего мне себя бояться? Я сам Ктулху местного значения. Так всем и говорю: "Не будите меня, я не Герцен, я Ктулху!". Даже песню сочинил:
http://samlib.ru/m/masherow_e_l/ctulhusong.shtml

Измерительные приборы обычно вообще в целом числе выдают, на уровне АЦП, имею в виду. И даже при скромном числе разрядов. Но если там, где "разность двух близких измерений" это предельная вычислительная проблема, можно в single поверить, то уже для скромного обращения матриц размером в сотни и более начинаются проблемы.

[Евгений Машеров]

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

[ewert]

Евгений Машеров в сообщении #1045017 писал(а):

Именно для такой функции имеет место простое рекуррентное соотношение для трёх последовательных точек,

Сугубо теорекхтицски. А пракхтицски -- Вы чересчур уж высокого мнения о точности отсчётов.

Евгений Машеров в сообщении #1044957 писал(а):

то уже для скромного обращения матриц размером в сотни и более начинаются проблемы.

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

[ewert]

Там, кстати, есть проблема. А что, если частота на много-многих периодах слежения плавает?... -- Тогда все уже предложенные рекомендации поплывут.

Это уже отдельная тема для разговора.

[Евгений Машеров]

Ну, "тут зависит". У задачи множество решений, и какое лучше - зависит от деталей, топикстартером не упомянутых.
Скажем, для "быстро и грубо" можно ограничится подсчётом переходов через $0$. Что даст точность определения частоты с погрешностью порядка процента или менее, а декремент затухания оценить, проведя огибающую через максимумы.
Если нужно повыжимать точность - нелинейная регрессия.
Если частота плавает - преобразование Гильберта. И затем
$f(t)+i\mathfrak{H}(f(t))=A(t)e^{i\varphi(t)}$
Затем находя переменную во времени амплитуду и мгновенную фазу, а от неё идя к частоте.
Что до того, что матрица $4$ на $4$, а не $100$ на $100$, то, увы, близости к вырожденной это не исключает, причём именно в этом случае производные будут коррелированы и дадут мультиколлинеарность со всеми вычислительными прелестями, так что я бы всё же пользовал double, оставив single для случая, когда надо хранить много, а считать мало.

[Mavlik]

Пока что не разобрался с методом. Нашел статью
http://www.hse.ru/data/2012/08/10/1256198474/Andreev_Applied%20Informatics_2010.pdf
с описанием алгоритма. Но как-то мне не очень понравилось про линеаризацию вектора-невязки.
Хотя, может быть, так и должно быть.
$y$ - вектор-столбец экспериментальных данных,
$f$ - функция
$r$ - столбец невязки
$c$ — вектор неизвестных параметров

$f - y=r$

$A^{(k)}\Delta c^{(k)}-(J^{T})^{(k)}r^{(k+1)}=-\upsilon ^{k}$
где $A^{(k)}=(J^{T})^{(k)}J^{(k)}$

В итоге получается вот такое выражение.
$(A^{(k)}+\lambda^{(k)}D)\Delta c^{(k)}=-\upsilon ^{(k)}$
$D=\operatorname{diag}(A^{(k)})$

Если верно понимаю, из него обращением матрицы $(A^{(k)}+\lambda^{(k)}D)$ находим $\Delta c^{(k)}$, а затем $c^{(k+1)}$
$c^{(k+1)}=c^{(k)}+\Delta c^{(k)}$

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

[Mavlik]

Но как-то не совсем понятно. Матрица Якоби непрямоугольная ? Тогда как взять якобиан?

[Евгений Машеров]

Зачем якобиан-то?