Вывод модели 2-гауссовой смеси с MCMC и PyMC

10

Проблема

Я хочу соответствовать модельным параметрам простой 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()

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

NLS подходят результаты

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

Байесовский подход: 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)

α

MCMC альфа-резюме

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

MCMC центров_0 резюме

α

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

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

user2304916
источник

Ответы:

6

Проблема вызвана тем, что 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'который дает плохие результаты (я не знаю, почему это так, но это, конечно, не похоже на разумное предложение для двоичной переменной).

Том Минка
источник
Спасибо, это действительно полезно, хотя это становится почти невыносимо медленным. Вы можете кратко объяснить смысл proposal_distributionи proposal_sdи почему использование Priorлучше для категориальных переменных?
user2304916
Спасибо, Том. Я согласен, что Пуассон здесь странный выбор. Я открыл вопрос: github.com/pymc-devs/pymc/issues/627
twiecki
2

σ

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 Книги ошибок они обсуждают этот тип модели и утверждают, что существуют проблемы сходимости без очевидного решения. Проверьте также здесь .

jpneto
источник
Это хороший момент, сейчас я использую гамму. Тем не менее, альфа-трасса всегда стремится к 0 (вместо 0,4). Интересно, есть ли какая-то глупая ошибка, скрывающаяся в моем примере.
user2304916
Я попробовал Гамму (.1, .1), но Exp (.1), кажется, работает лучше. Кроме того, когда автокорреляция высока, вы можете включить некоторое прореживание, например,mcmc.sample(60000, 3000, 3)
jpneto
2

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

Бенджамин Даути
источник