2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки


Правила форума


В этом разделе нельзя создавать новые темы.

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

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

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

Обязательно просмотрите тему Правила данного раздела, иначе Ваша тема может быть удалена или перемещена в Карантин, а Вы так и не узнаете, почему.



Начать новую тему Ответить на тему На страницу Пред.  1, 2, 3  След.
 
 Re: Подгонка функции (curve fitting)
Сообщение11.08.2015, 19:38 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение12.08.2015, 12:15 
Заслуженный участник


11/05/08
32166
Mavlik в сообщении #1044461 писал(а):
$\alpha=\frac{1}{_{t_{1}}}\cdot \ln(\frac{sz}{\sin (\omega\cdot t_{1})})$

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

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

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

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

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

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

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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение12.08.2015, 21:09 
Супермодератор
Аватара пользователя


20/11/12
5728
 ! 
Евгений Машеров в сообщении #1044570 писал(а):
Матрица одна. nxn.
Евгений Машеров, замечание за хроническое неоформление формул. В следующий раз буду сносить ваши посты в Карантин.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение12.08.2015, 21:53 
Заслуженный участник
Аватара пользователя


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

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

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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение12.08.2015, 22:56 
Заслуженный участник


11/05/08
32166
Евгений Машеров в сообщении #1044858 писал(а):
А здесь лишь указание на размерность.

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

Евгений Машеров в сообщении #1044858 писал(а):
Затрудняюсь сказать, где он может быть полезен - для хранения больших массивов при острой нехватке памяти и отсутствии больших вычислений с ними?

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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение13.08.2015, 09:18 
Заслуженный участник
Аватара пользователя


11/03/08
9490
Москва

(Оффтоп)

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


Давненько я не брал в руки 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$
(насколько малые - ну, вот я увеличивал диагональные на один процент и сверх того прибавлял одну миллионную).
Скорректировав значения параметров - повторяем расчёт, следя за сходимостью.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение13.08.2015, 11:27 
Заслуженный участник
Аватара пользователя


01/08/06
3049
Уфа
Евгений Машеров в сообщении #1044858 писал(а):
Тип single, 32-битный плавающий, вообще стоит избегать. Затрудняюсь сказать, где он может быть полезен
Побойтесь ктулху! Он отлично подходит почти для всех исходных данных. Редко где измерительные приборы могут обеспечить точность больше 7 значащих цифр. Да и простые операции с исходными данными (2-3 штуки) чаще всего не приводят к неприемлемому уменьшению точности. Я с геофизическими данными работаю, где почти всё в single. И даже там, где нужно получить разность двух близких измерений (частая операция), single всё ещё отлично работает.

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение13.08.2015, 11:50 
Заслуженный участник
Аватара пользователя


11/03/08
9490
Москва
worm2 в сообщении #1044946 писал(а):
Евгений Машеров в сообщении #1044858 писал(а):
Тип single, 32-битный плавающий, вообще стоит избегать. Затрудняюсь сказать, где он может быть полезен
Побойтесь ктулху! Он отлично подходит почти для всех исходных данных. Редко где измерительные приборы могут обеспечить точность больше 7 значащих цифр. Да и простые операции с исходными данными (2-3 штуки) чаще всего не приводят к неприемлемому уменьшению точности. Я с геофизическими данными работаю, где почти всё в single. И даже там, где нужно получить разность двух близких измерений (частая операция), single всё ещё отлично работает.

(Оффтоп)

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


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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение13.08.2015, 14:53 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение13.08.2015, 22:36 
Заслуженный участник


11/05/08
32166
Евгений Машеров в сообщении #1045017 писал(а):
Именно для такой функции имеет место простое рекуррентное соотношение для трёх последовательных точек,

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

Евгений Машеров в сообщении #1044957 писал(а):
то уже для скромного обращения матриц размером в сотни и более начинаются проблемы.

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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение14.08.2015, 00:57 
Заслуженный участник


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

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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение14.08.2015, 09:36 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение17.08.2015, 11:28 


27/03/12
23
Пока что не разобрался с методом. Нашел статью
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)}$

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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение19.08.2015, 08:47 


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

 Профиль  
                  
 
 Re: Подгонка функции (curve fitting)
Сообщение19.08.2015, 10:12 
Заслуженный участник
Аватара пользователя


11/03/08
9490
Москва
Зачем якобиан-то?

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 32 ]  На страницу Пред.  1, 2, 3  След.

Модераторы: Модераторы Математики, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group