Проблема

Я хочу соответствовать модельным параметрам простой 2-гауссовой смеси населения. Учитывая всю шумиху вокруг байесовских методов, я хочу понять, является ли для этой проблемы байесовский вывод лучшим инструментом, чем традиционные методы подбора.

Пока MCMC работает очень плохо в этом игрушечном примере, но, возможно, я просто что-то упустил. Итак, давайте посмотрим код.

Инструменты

Я буду использовать python (2.7) + стек scipy, lmfit 0.8 и PyMC 2.3.

Блокнот для воспроизведения анализа можно найти здесь

Генерация данных

Сначала давайте сгенерируем данные:

from scipy.stats import distributions

# Sample parameters
nsamples = 1000
mu1_true = 0.3
mu2_true = 0.55
sig1_true = 0.08
sig2_true = 0.12
a_true = 0.4

# Samples generation
np.random.seed(3)  # for repeatability
s1 = distributions.norm.rvs(mu1_true, sig1_true, size=round(a_true*nsamples))
s2 = distributions.norm.rvs(mu2_true, sig2_true, size=round((1-a_true)*nsamples))
samples = np.hstack([s1, s2])

Гистограмма samplesвыглядит так:

«широкий пик», компоненты трудно определить на глаз.

Классический подход: подгонка гистограммы

Давайте сначала попробуем классический подход. Используя lmfit , легко определить 2- пиковую модель:

import lmfit

peak1 = lmfit.models.GaussianModel(prefix='p1_')
peak2 = lmfit.models.GaussianModel(prefix='p2_')
model = peak1 + peak2

model.set_param_hint('p1_center', value=0.2, min=-1, max=2)
model.set_param_hint('p2_center', value=0.5, min=-1, max=2)
model.set_param_hint('p1_sigma', value=0.1, min=0.01, max=0.3)
model.set_param_hint('p2_sigma', value=0.1, min=0.01, max=0.3)
model.set_param_hint('p1_amplitude', value=1, min=0.0, max=1)
model.set_param_hint('p2_amplitude', expr='1 - p1_amplitude')
name = '2-gaussians'

Наконец мы подгоняем модель с помощью симплексного алгоритма:

fit_res = model.fit(data, x=x_data, method='nelder')
print fit_res.fit_report()

В результате получается следующее изображение (красные пунктирные линии соответствуют центрам):

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

Байесовский подход: MCMC

Я определяю модель в PyMC иерархически. centersи sigmasраспределение априоров для гиперпараметров, представляющих 2 центра и 2 сигмы 2 гауссианов. alphaэто доля первой популяции, а предыдущее распределение здесь бета.

Категориальная переменная выбирает между двумя популяциями. Насколько я понимаю, эта переменная должна быть того же размера, что и data ( samples).

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

sigmas = pm.Normal('sigmas', mu=0.1, tau=1000, size=2)
centers = pm.Normal('centers', [0.3, 0.7], [1/(0.1)**2, 1/(0.1)**2], size=2)
#centers = pm.Uniform('centers', 0, 1, size=2)

alpha  = pm.Beta('alpha', alpha=2, beta=3)
category = pm.Categorical("category", [alpha, 1 - alpha], size=nsamples)

@pm.deterministic
def mu(category=category, centers=centers):
    return centers[category]

@pm.deterministic
def tau(category=category, sigmas=sigmas):
    return 1/(sigmas[category]**2)

observations = pm.Normal('samples_model', mu=mu, tau=tau, value=samples, observed=True)
model = pm.Model([observations, mu, tau, category, alpha, sigmas, centers])

Затем я запускаю MCMC с довольно большим количеством итераций (1e5, ~ 60 с на моей машине):

mcmc = pm.MCMC(model)
mcmc.sample(100000, 30000)

$\alpha$

Также центры гауссианов также не сходятся. Например:

$\alpha$

Так что здесь происходит? Я делаю что-то не так или MCMC не подходит для этой проблемы?

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

Проблема вызвана тем, что PyMC рисует образцы для этой модели. Как объяснено в разделе 5.8.1 документации PyMC, все элементы переменной массива обновляются вместе. Для таких маленьких массивов centerэто не проблема, но для большого массива categoryэто приводит к низкой скорости приема. Вы можете увидеть уровень принятия через

print mcmc.step_method_dict[category][0].ratio

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

import pymc as pm
sigmas = pm.Normal('sigmas', mu=0.1, tau=1000, size=2)
centers = pm.Normal('centers', [0.3, 0.7], [1/(0.1)**2, 1/(0.1)**2], size=2)
alpha  = pm.Beta('alpha', alpha=2, beta=3)
category = pm.Container([pm.Categorical("category%i" % i, [alpha, 1 - alpha]) 
                         for i in range(nsamples)])
observations = pm.Container([pm.Normal('samples_model%i' % i, 
                   mu=centers[category[i]], tau=1/(sigmas[category[i]]**2), 
                   value=samples[i], observed=True) for i in range(nsamples)])
model = pm.Model([observations, category, alpha, sigmas, centers])
mcmc = pm.MCMC(model)
# initialize in a good place to reduce the number of steps required
centers.value = [mu1_true, mu2_true]
# set a custom proposal for centers, since the default is bad
mcmc.use_step_method(pm.Metropolis, centers, proposal_sd=sig1_true/np.sqrt(nsamples))
# set a custom proposal for category, since the default is bad
for i in range(nsamples):
    mcmc.use_step_method(pm.DiscreteMetropolis, category[i], proposal_distribution='Prior')
mcmc.sample(100)  # beware sampling takes much longer now
# check the acceptance rates
print mcmc.step_method_dict[category[0]][0].ratio
print mcmc.step_method_dict[centers][0].ratio
print mcmc.step_method_dict[alpha][0].ratio

proposal_sdИ proposal_distributionопции описаны в разделе 5.7.1 . Для центров я установил предложение примерно равным стандартному отклонению апостериорного значения, которое намного меньше значения по умолчанию из-за объема данных. PyMC пытается настроить ширину предложения, но это работает только в том случае, если ваш уровень принятия достаточно высок для начала. Ибо category, по умолчанию, proposal_distribution = 'Poisson'который дает плохие результаты (я не знаю, почему это так, но это, конечно, не похоже на разумное предложение для двоичной переменной).

$\sigma$

sigmas = pm.Exponential('sigmas', 0.1, size=2)

Обновить:

Я получил около начальных параметров данных, изменив эти части вашей модели:

sigmas = pm.Exponential('sigmas', 0.1, size=2)
alpha  = pm.Beta('alpha', alpha=1, beta=1)

и вызывая mcmc с некоторым прореживанием:

mcmc.sample(200000, 3000, 10)

Результаты:

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

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