Установка бимодального гауссовского распределения с некоторыми фиксированными параметрами

Question

Задача: Я хочу сопоставить эмпирические данные с бимодальным нормальным распределением, из которого я знаю из физического контекста расстояние от пиков (фиксированное), а также то, что оба пика должны иметь такое же стандартное отклонение.

Я пытался создать собственный дистрибутив с scipy.stats.rv_continous (см. Код ниже), но параметры всегда соответствуют 1. Кто-то понимает, что происходит, или может указать мне другой подход к решению проблемы?

Детали: Я избегал loc и scale параметров и реализовать их как m и s непосредственно в _pdf -метод, так как пиковое расстояние delta не должно зависеть от scale. Чтобы компенсировать это, я установил их floc=0 и fscale=1 в fit -методе и на самом деле хочу FIT-параметры для m, s и весы вершин w

Что ожидать в данной выборке является распределением с пиками вокруг x=-450 и x=450 (=>m=0). Stdev s должен составлять около 100 или 200, но не 1,0, а вес w должен составлять ок. 0,5

from __future__ import division 
from scipy.stats import rv_continuous 
import numpy as np 


class norm2_gen(rv_continuous): 
    def _argcheck(self, *args): 
     return True 

    def _pdf(self, x, m, s, w, delta): 
     return np.exp(-(x-m+delta/2)**2/(2. * s**2))/np.sqrt(2. * np.pi * s**2) * w + \ 
       np.exp(-(x-m-delta/2)**2/(2. * s**2))/np.sqrt(2. * np.pi * s**2) * (1 - w) 


norm2 = norm2_gen(name='norm2') 

data = [487.0, -325.5, -159.0, 326.5, 538.0, 552.0, 563.0, -156.0, 545.5, 341.0, 530.0, -156.0, 473.0, 328.0, -319.5, -287.0, -294.5, 153.5, -512.0, 386.0, -129.0, -432.5, -382.0, -346.5, 349.0, 391.0, 299.0, 364.0, -283.0, 562.5, -42.0, 214.0, -389.0, 42.5, 259.5, -302.5, 330.5, -338.0, 508.5, 319.5, -356.5, 421.5, 543.0] 

m, s, w, delta, loc, scale = norm2.fit(data, fdelta=900, floc=0, fscale=1) 
print m, s, w, delta, loc, scale 

>>> 1.0 1.0 1.0 900 0 1 

Answer 1

Я был в состоянии получить распределение по размеру данных после того, как сделать пару настроек:

• Использование w, как вы делали, у вас есть подразумеваемое ограничение, 0 < = w < = 1 Решатель, используемый методом fit(), не знает об этом ограничении, поэтому w может привести к необоснованным значениям. Один из способов обращения с этим типом ограничений состоит в том, чтобы позволить w быть произвольным реальным значением, но в формуле для PDF конвертировать w в долю phi между 0 и 1 с использованием phi = 0.5 + arctan(w)/pi.
• В методе общего назначения fit() используется процедура численной оптимизации для определения оценки максимального правдоподобия. Как и большинство таких подпрограмм, для его оптимизации нужна начальная точка. Начальная точка по умолчанию - все 1 с, но это не всегда хорошо работает. Вы можете выбрать другую отправную точку, предоставив значения в качестве позиционных аргументов fit() после данных. Значения, которые я использовал в скрипте, работали; Я не изучил, насколько чувствительны результаты этих стартовых значений.

Я сделал две оценки. Во-первых, я позволил delta быть свободным параметром, а во втором, я установил delta на 900.

Приведенный ниже сценарий генерирует следующий сюжет:

Вот сценарий:

from __future__ import division 
from scipy.stats import rv_continuous 
import numpy as np 
import matplotlib.pyplot as plt 


class norm2_gen(rv_continuous): 
    def _argcheck(self, *args): 
     return True 

    def _pdf(self, x, m, s, w, delta): 
     phi = 0.5 + np.arctan(w)/np.pi 
     return np.exp(-(x-m+delta/2)**2/(2. * s**2))/np.sqrt(2. * np.pi * s**2) * phi + \ 
       np.exp(-(x-m-delta/2)**2/(2. * s**2))/np.sqrt(2. * np.pi * s**2) * (1 - phi) 

norm2 = norm2_gen(name='norm2') 


data = [487.0, -325.5, -159.0, 326.5, 538.0, 552.0, 563.0, -156.0, 545.5, 
     341.0, 530.0, -156.0, 473.0, 328.0, -319.5, -287.0, -294.5, 153.5, 
     -512.0, 386.0, -129.0, -432.5, -382.0, -346.5, 349.0, 391.0, 299.0, 
     364.0, -283.0, 562.5, -42.0, 214.0, -389.0, 42.5, 259.5, -302.5, 
     330.5, -338.0, 508.5, 319.5, -356.5, 421.5, 543.0] 

# In the fit method, the positional arguments after data are the initial 
# guesses that are passed to the optimization routine that computes the MLE. 
# First let's see what we get if delta is not fixed. 
m, s, w, delta, loc, scale = norm2.fit(data, 1.0, 1.0, 0.0, 900.0, floc=0, fscale=1) 

# Fit the disribution with delta fixed. 
fdelta = 900 
m1, s1, w1, delta1, loc, scale = norm2.fit(data, 1.0, 1.0, 0.0, fdelta=fdelta, floc=0, fscale=1) 

plt.hist(data, bins=12, normed=True, color='c', alpha=0.65) 
q = np.linspace(-800, 800, 1000) 
p = norm2.pdf(q, m, s, w, delta) 
p1 = norm2.pdf(q, m1, s1, w1, fdelta) 
plt.plot(q, p, 'k', linewidth=2.5, label='delta=%6.2f (fit)' % delta) 
plt.plot(q, p1, 'k--', linewidth=2.5, label='delta=%6.2f (fixed)' % fdelta) 
plt.legend(loc='best') 
plt.show()