Как определить параметры синуса?

Есть результаты физического эксперимента, которые с некоторой ошибкой могут быть представлены в виде функции
y = A*sin(B*x + x0) + y0
Всего пар x;y около 1000, причем в эксперимент попадает несколько(8-10) периодов функции.
Иными словами точки насажены довольно плотно.

Вопрос - каким образом можно определить A, B, x0, y0 так, что бы точность определения параметров росла с ростом числа известных точек, ну и желательно что бы минимизировалась ошибка(но уже необязательно)?

Это называется "нелинейный МНК". Используй Левенберга-Марквардта для минимизации суммарной ошибки, т.е. функции вида F=(y(x0)-y0)^2 + ... + (y(xn)-yn)^2

Реализация метода - либо пакет levmar (c/c++), либо из ALGLIB - http://alglib.sources.ru/optimization/levenbergmarquardt.php

Добавлено 11.10.09, 17:04
Добавление 1: да, минимизация ведется по A,B, xo, yo.

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

(спасибо, пока думаю)

Добавлено 11.10.09, 19:13
Тут нужна некоторая помощь.

У нас выборка x(i), y(i)
Как я понимаю, в качестве f(i) берется (Y - A*sin(B*x(i) + X0) + Y0)^2
Где Y,A,X0,Y0 - переменные от которых зависит функция.
Я взял реализацию на C# и применял ее вот так:

        Dim s As New alglib.minlm.lmstate
        Dim sol() As Double = {0.25, 1591.549431, 0, -150} ' {0.1, 1000, 0, -120} '  'Примерно подобранные начальные значения
        'ReDim sol(results.Count)
        alglib.minlm.minlmfj(4, results.Count, sol, 0.0000001, 0.0000001, 10000, s)' Тут числа взяты наугад
 
        While alglib.minlm.minlmiteration(s)
            Dim A As Double = s.x(0)
            Dim B As Double = s.x(1)
            Dim X0 As Double = s.x(2)
            Dim Y0 As Double = s.x(3)
            If s.needf Then
                Dim _F As Double = 0
                For Each R As PointF In results
                    _F += (R.Y - (A * Math.Sin(B * R.X + X0) + Y0)) ^ 2
                Next
                s.f = _F
                Debug.Print(_F.ToString) 'ТУТ вывожу значения вычисляемой функции
            End If
            If s.needfij Then
                For i As Integer = 0 To results.Count - 1
                    Dim r As PointF = results(i)
                    s.fi(i) = r.Y - (A * Math.Sin(B * r.X + X0) + Y0)
                    s.j(i, 0) = Math.Sin(B * r.X + X0) REM dF/dA
                    s.j(i, 1) = A * B * Math.Cos(B * r.X + X0) REM dF/dB
                    s.j(i, 2) = A * Math.Cos(B * r.X + X0) REM dF/dX0
                    s.j(i, 3) = 1 REM dF/dY0
                Next
            End If
            If s.needfg Or s.needfgh Then
                Throw New Exception()
            End If
        End While
        Dim rep As alglib.minlm.lmreport
        Dim X(4) As Double
        alglib.minlm.minlmresults(s, X, rep)
        Dim z As String = "{" & X(0) & ", " & X(1) & ", " & X(2) & ", " & X(3) & "}"
        MsgBox(z)

Использую вот такой код.


Теперь пробую прогнать ее на одной и той же выборке, при разных значениях стартовых параметров
При стартовых параметрах {0.25, 1591.549431, 0, -150} (это довольно близко к правде) алгоритм дает ответ
{0.250000043588598, 1591.54946567717, 2.17883057160585E-08, -150.00006958343}
Что похоже на правду.
При этом последовательные значения функции вот такие:

Цитата

175285943.807378
175279093.429274
174998443.565119
155484743.497379
47898495.6142789
43850845.497696

Если же взять другие начальные значения, например {0.1, 1000, 0, -120} что уже довольно далеко от истины, получается ответ
{0.0999999859957626, 1000.00000574627, 5.74626982824354E-09, -120.000062631127}
А это совсем неверно
При этом вычисленные значения функции следующие:

Цитата

142013745.650736
142008195.060041
141780793.005068
125971771.737748
38809425.2111686
35530193.4296744

Вообще какие-то громадные числа получаются. Причина завершения в обоих случаях "relative step is no more than EpsX"

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

Гм. Ну а как тогда все таки решать исходную задачу? Нужно найти, получается, глобальный минимум, причем в нем f(x) = 0
И посмотрите плз, правильно ли посчитан якобиан

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

Что?)

Э-э-э... может я не понял суть задачи .. в более понятных мне терминах радиотехнических, она звучит так: необходимо определить амплитуду, частоту, начальную фазу и постоянную составляющую синусоидального сигнала, заданного дискретными значениями. Только значения y должны быть "взяты" через равные промежутки времени, а в Вашей формуле вместо времени фигурирует некий x, и о нем ничего неизвестно. Если растояния между соседними значениями x одинаковы, можно вычислить спектр сигнала.

Сообщение отредактировано: Prince - 12.10.09, 10:20

Цитата

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

Да, совершенно верно, я о том же подумал.
Левенберг-Марквардт в данном случае - это чересчур.

ANDLL, неправильно считается Якобиан.

Цитата ANDLL @ 11.10.09, 18:49

s.j(i, 0) = Math.Sin(B * r.X + X0) REM dF/dA
s.j(i, 1) = A * B * Math.Cos(B * r.X + X0) REM dF/dB
s.j(i, 2) = A * Math.Cos(B * r.X + X0) REM dF/dX0
s.j(i, 3) = 1 REM dF/dY0

Во-первых, в каждой строчке потерян минус. Во-вторых, в строке "s.j(i,1)=" есть лишний множитель B.

Гм.
Только наверное
s.j(i, 1) = - A * r.X * Math.Cos(B * r.X + X0) REM dF/dB

Цитата 0797293 @ 12.10.09, 10:21

Цитата
А нельзя ли для нахождения параметров синуса применить спектральный анализ?

Да, совершенно верно, я о том же подумал.
Левенберг-Марквардт в данном случае - это чересчур.

Мне кажется, это как раз спектральный анализ - не в тему Во-первых, методы, основанные на БПФ, позволяют работать только с точками на равномерной сетке - а здесь сетка неравномерная. Во-вторых, спектральные методы позволяют выбирать множитель "B" только из конечного диапазона частот, а Левенберг-Марквардт позволяет непрерывно выбирать оптимальный параметр. В-третьих, спектральный анализ совершенно не годится для ситуации, когда нельзя однозначно определить значение функции в точке (есть несколько точек с совпадающими абсциссами, но разными значениями функции) - ситуация, часто встречающаяся в практике.


P.S. Вообще, наверное надо добавить в алгоритм какой-то код для проверки точности вычисления Якобиана/градиента/Гессиана.

Добавлено 12.10.09, 10:28

Цитата ANDLL @ 12.10.09, 10:23

Гм.
Только наверное
s.j(i, 1) = - A * r.X * Math.Cos(B * r.X + X0) REM dF/dB

А, ну да, конечно же вы правы я неясно выразился.

Заменил код.
Вот так правильный якобиан?

                    s.j(i, 0) = -Math.Sin(B * r.X + X0) REM dF/dA
                    s.j(i, 1) = -A * r.X * Math.Cos(B * r.X + X0) REM dF/dB
                    s.j(i, 2) = -A * Math.Cos(B * r.X + X0) REM dF/dX0
                    s.j(i, 3) = -1 REM dF/dY0

А то результат все равно довольно неправильный(

Цитата

Во-первых, методы, основанные на БПФ, позволяют работать только с точками на равномерной сетке

Вот в этом и была неясность. Раз сетка неравномерная, со спектром я промахнулся.

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