Gaussian fit

[Adventor]

Здравствуйте!

Есть такая задача: дан массив. Нужно аппроксимировать его с помошью гауссиан. Известно количество гауссиан и положение их центров.


Пример результата:


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

Для данного примера я вижу что следует минимизировать функционал зависящий от 6 параметров:
$y (A_1, \sigma_1, A_2, \sigma_2, A_3, \sigma_3 ) = A_1 e^{ {-}\frac{(x - \mu_1)^2}{(2 \sigma_1)^2} } + A_2 e^{ {-} \frac{(x - \mu_2)^2}{(2 \sigma_2)^2} } + A_3 e^{ {-} \frac{(x - \mu_3)^2}{(2 \sigma_3)^2} }$

Наивный вариант: задать сетку значений всех 6 параметров, и найти минимум для всех вариантов параметров. Но встает вопрос как выбирать оптимальный шаг на сетке параметров и какая будет вычислительная сложность. Стоит ли так делать?

[Pphantom]

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

А так... метод Левенберга-Марквардта не годится?

[Adventor]

Pphantom в сообщении #1330473 писал(а):

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

А так... метод Левенберга-Марквардта не годится?

оупс, конечно.

Надо подумать...

[dsge]

В функционале должна быть сумма по всем наблюдениям. Сумма $A_i$ должна равняться единице, и они должны быть положительными, иначе это не будет плотность. Т.е. оптимизация условная. Далее любой солвер-оптимизатор найдет максимум правдоподобия. Функционал надо максимизировать. Чтобы был «честный» максимум правдоподобия, экспоненты надо разделить на сигмы и корень из два пи.

[Adventor]

dsge в сообщении #1330500 писал(а):

В функционале должна быть сумма по всем наблюдениям. Сумма $A_i$ должна равняться единице, и они должны быть положительными, иначе это не будет плотность. Т.е. оптимизация условная. Далее любой солвер-оптимизатор найдет максимум правдоподобия. Функционал надо максимизировать. Чтобы был «честный» максимум правдоподобия, экспоненты надо разделить на сигмы и корень из два пи.

В моем случае гауссианы не нормализованы на единицу.

[dsge]

Надо не плотность подогнать, а данные? Тогда, да.

[Adventor]

dsge в сообщении #1330584 писал(а):

Надо не плотность подогнать, а данные? Тогда, да.

Плотность вообще никак не фигурирует.

[photon]

Гуглите GMM algorithm и обрящете искомое. (GMM = Gaussian Mixture Model). У вас, правда, задача попроще - количество и центры уже известны, но упрощать - не усложнять.

[B@R5uk]

Adventor в сообщении #1330471 писал(а):

Для данного примера я вижу что следует минимизировать функционал зависящий от 6 параметров:

$$y (A_1, \sigma_1, A_2, \sigma_2, A_3, \sigma_3 ) = A_1 e^{ {-}\frac{(x - \mu_1)^2}{(2 \sigma_1)^2} } + A_2 e^{ {-} \frac{(x - \mu_2)^2}{(2 \sigma_2)^2} } + A_3 e^{ {-} \frac{(x - \mu_3)^2}{(2 \sigma_3)^2} }$$

По параметрам $A_i$ оптимизацию можно сделать сразу аналитически. То есть численная оптимизация будет проходить только по трём параметрам: $\sigma_i$.

Если вас интересует простой надёжный алгоритм, не использующий производные, то рекомендую алгоритм Нелдера-Мида. Именно он используется в стандартной функции fminsearch программы МАТЛАБ. Его простоту гарантирую, так как сам его реализовывал на си.

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

Чтобы лучше понять в чём именно заключается проблема, представьте себе, что вам необходимо разложить трёхмерный вектор по неортогональному базису в трёхмерном пространстве. Причём два вектора 1-й и 2-й этого базиса практически параллельны. Тогда небольшая флуктуация координат раскладываемого вектора приведёт к гигантской флуктуации коэффициентов разложения по 1-ому и 2-ому вектору. Аналогичная ситуация возникает и у вас, так как направления в пространстве допустимых моделей, задаваемые первой и второй гауссианиной практически параллельны.

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