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

[Mavlik]

Есть такой вопрос.
У меня имеется осциллограмма(набор точек в файле) - сигнал с резистора в цепи RLC при разрядке конденсатора.
Сигнал представляет затухающую синусоиду.
$y(t)=I_{0}\cdot e^{-\alpha \cdot t}\sin (\omega\cdot t)$
Неизвестные коэффициенты:
$I_{0}, \alpha, \omega$
Точек много, но для решения задачи можно брать не все. а, скажем, каждую сотую точку.
Пробовал решать с помощью МНК, но застрял после того, как продифференцировал функцию невязок.
(Как ее дальше решать, я не знаю. А тем более программно).

$\sigma(I_{0},\alpha,\omega)=\sum_{i=1}^{n}(y_{i}-I_{0}\cdot e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))^{2}$
Получается вот такая система уравнений.
$\frac{\partial}{\partial I_{0}}\sigma(I_{0},\alpha,\omega)=-2\sum_{i=1}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))\cdot e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i})=0$
$\frac{\partial}{\partial \alpha}\sigma(I_{0},\alpha,\omega)=2\sum_{i=1}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))\cdot I_{0}\cdot t_{i}\cdot e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i})=0$
$\frac{\partial}{\partial \omega}\sigma(I_{0},\alpha,\omega)=-2\sum_{i=1}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))\cdot I_{0}\cdot t_{i}\cdot e^{-\alpha \cdot t_{i}}\cos (\omega\cdot t_{i})=0$

Я хочу написать программу на языке программирования, а не просто решить в каком-то пакете.
В одном из математических пакетов есть NonLinear Curve Fit, там делается большое число итераций, после чего считаются коэффициенты. Вроде бы используется метод Левенберга-Марквардта.
Буду благодарен за подсказку в том, как решить поставленную выше задачу.

[Александрович]

Выпрямить и по максимальным и минимальным точкам все найти.

[Mavlik]

что значит "выпрямить"?
Нужно именно для данного набора точек.

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

Если точек достаточно много (ну, скажем, с десяток на период), то максимум (минимум), определяемый, как $x_{i-1}\le x_i\ge x_{i+1}$ и аналогично для минимума, будет достаточно близок к максимальным/минимальным значениям функции в окрестности точки i. Тогда, составив последовательность из максимумов и (абсолютных величин) минимумов и соответствующих им моментов времени, и прологарифмировав значения минимумов и максимумов, получим линейную регрессию, в которой коэффициент при времени будет соответствовать показателю экспоненты. А разделив период времени на число минимумов и максимумов в нём, получим период колебаний (половину периода, уточняет сержант-вычислитель Равно из батареи капитана Очевилность). Другая оценка периода - через число пересечений нуля.
Это довольно грубая оценка, но может оказаться достаточно точной для практики. Грубость её связана прежде всего с тем, что отобранные так точки выборки не попадают в точности на максимумы и минимумы, хотя, если их много, лежат достаточно близко и численно близки к максимумам и минимумам функции. Кроме того, принимается, что в максимуме синус, как сомножитель функции, равен единице (минус единице в минимуме), но ввиду экспоненциального множителя это не так, хотя при небольшом показателе экспоненты это работает, как приближение.
Даже если точность такой оценки недостаточна, она может быть использована, как начальное приближение (в том же Левенберге-Марквардте)

-- 07 авг 2015, 12:17 --

Теперь по Левенбергу-Марквардту.
Начинаем с какого-то приближения, получаем оценку $y(t)$, подставив в выражение приближенные значения параметров, находим отклонения фактических значений от полученных.
Затем пользуемся линейным приближением для этого отклонения
$\Delta y=Z\Delta a+\varepsilon$
где $\Delta a$ - отклонения истинных значений параметров от приближённых, $\varepsilon$ - включает как собственно ошибки измерения и пр., так и влияние неточности линейной модели, Z - матрица из значений производных в точках.
Далее обычная линейная регрессия даёт $\Delta a=(Z^TZ)^{-1}Z^T\Delta y$, позволяя уточнить значение вектора параметров, но значения производных обычно коррелированы, и обращать придётся почти вырожденную матрицу. Поэтому к её диагонали прибавляют константу $\lambda$ (метод Левенберга) или умножают диагональный элемент на $(1+\lambda)$ (собственно метод Левенберга-Марквардта). Чем больше лямбда - тем устойчивее численно, но медленнее (при очень большой получается метод наискорейшего спуска)

[Mavlik]

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

$I_{0}^{(k+1)}=\varphi_{1}(I_{0}^{(k)}, \alpha^{(k)}, \omega^{(k)})$
$\alpha^{(k+1)}=\varphi_{2}(I_{0}^{(k)}, \alpha^{(k)}, \omega^{(k)})$
$\omega^{(k+1)}=\varphi_{3}(I_{0}^{(k)}, \alpha^{(k)}, \omega^{(k)})$

[Алексей К.]

Мои правки к тому, что Mavlik в сообщении #1043204 писал(а):

Сигнал представляет затухающую синусоиду.
$y(t)=I_{0}\cdot e^{-\alpha \cdot t}\sin (\omega\, t{\color{magenta}+\varphi})$
Неизвестные коэффициенты:
$I_{0}, \alpha, \omega, \color{magenta}\varphi$

По-моему, нельзя для таких данных игнорировать сдвиг, 4-й параметр модели.

[Mavlik]

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

[Александрович]

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

что значит "выпрямить"?
Нужно именно для данного набора точек.

Выпрямить, значит взять модуль.

[Александрович]

Алексей К. в сообщении #1043277 писал(а):

Мои правки к тому, что Mavlik в сообщении #1043204 писал(а):

Сигнал представляет затухающую синусоиду.
$y(t)=I_{0}\cdot e^{-\alpha \cdot t}\sin (\omega\, t{\color{magenta}+\varphi})$
Неизвестные коэффициенты:
$I_{0}, \alpha, \omega, \color{magenta}\varphi$

По-моему, нельзя для таких данных игнорировать сдвиг, 4-й параметр модели.

Можно начать с первого нулевого значения функции. А сдвиг по времени учесть потом при нахождении $I_0$.

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

Со сдвигом - я бы рекомендовал перепараметризовать, используя тождество $A\sin(\omega t+\varphi)=B_0\sin \omega t+ B_1\cos \omega t$

[Mavlik]

Можно из вышеприведенной системы выразить первое слагаемое суммы и из него взять выражение для каждой переменной. Обозначим некоторый повторяющийся множитель как
$sz$

Далее вот так

$sz=\frac{\sum_{i=2}^{n}(y_{i}-I_{0}e^{-\alpha \cdot t_{i}}\sin (\omega\cdot t_{i}))}{I_{0}e^{-\alpha \cdot t_{1}}\sin (\omega\cdot t_{1})-y_{1}}$

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

$I_{0}=\frac{sz}{sin (\omega\cdot t_{1})\cdot t_{1}\cdot e^{-\alpha \cdot t_{1}}}$

$\omega=\frac{1}{t_{1}}\cdot \arccos\left ( \frac{sz}{sin (\omega\cdot t_{1})\cdot I_{0}\cdot t_{1}\cdot e^{-\alpha \cdot t_{1}}} \right )$

А вот программа на паскале. Но вылетает. Возможно, метод так нельзя применять.

program mnk;

uses
  SysUtils,math;

var
  xt,yt:array[1..5000] of single;
  x,dx:single;
  I0,omega,alpha:single;
  I0k,omegak,alphak:array[1..5000] of single;
  i,j,n:integer;
  sum,sum1:single;
  sz,znam:single;

begin
I0:=2.4;
alpha:=0.1;
omega:=3.3;
x:=0;
dx:=0.1;
n:=100;

for i:=1 to n do
begin
x:=x+dx;
xt[i]:=x;
yt[i]:=I0*exp(-alpha*xt[i])*sin(omega*xt[i]);
end;

I0k[1]:=10;
alphak[1]:=6.12;
omegak[1]:=1;

for j:=1 to 1000 do
begin
 sum1:=0;
 for i:=2 to n do
  begin
   sum1:=sum1+(yt[i]-I0k[j]*exp(-alphak[j]*xt[i])*sin(omegak[j]*xt[i]));
  end;
znam:=-(yt[1]-I0k[j]*exp(-alphak[j]*xt[1]));
sz:=sum1/znam;

alphak[j+1]:= (-1/xt[1])*ln(abs(sz/sin(omegak[j]*xt[1])));
I0k[j+1]:=sz/xt[1]/exp(-alphak[j]*xt[1]);//sin(omegak[j]*xt[1]);
omegak[j+1]:=(1/xt[1])*arccos(sz/I0k[j]/xt[1]/exp(-alphak[j]*xt[1]));

end;
end.

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

0. Я бы рекомендовал всё же включить в модель фазовый сдвиг, как ещё один параметр. Даже если Вы можете гарантировать, что в начальный момент аргумент синуса равен нулю - начальный момент процесса и первое снятие его значения не обязаны совпадать.
1. Я не вижу никаких оснований ожидать того, что выписанные Вами выражения будут сходиться. Метод простой итерации сходится при довольно узких условиях.
2. Скажите, а что будет, если Вы в выражение для sz подставите точные значения параметров, а ошибка при первом измерении будет 0? И что будет, если ошибка при нём не 0?
2. А что сложного в Марквардте? Производные Вы считать можете, перемножить матрицы недорого, обратить немногим сложнее.

[Mavlik]

Евгений Машеров
Я не знаю, что будет при подстановке точных значений.
А метод ньютона тоже не подойдет?

Где можно почитать про метод марквардта?

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

0. Если вдруг подставленные в формулу параметры окажутся равны точным значениям, а ошибка измерения при первом наблюдении окажется нулевой - у Вас в выражении для sz будет деление на ноль. А если ошибка не ноль - то деление на величину ошибки. То есть результат закинет неизвестно куда.
1. Ньютона - подойдёт. Только если Вы не будете учитывать корреляции между значениями производных, то сходиться будет медленно и печально. Или не будет. А если учтёте - получится Левенберг-Марквардт.
2. Лично я прочёл в книге Е.З.Демиденко "Линейная и нелинейная регрессии". Но он во многих местах описан.

[Mavlik]

Вот, нашел описание метода.
http://optimizaciya-sapr.narod.ru/bez_mnogomer/makva.html
Получается, что нужно решать каждое уравнение системы, и для каждого будет своя матрица Гессе по всем четырем переменным?