O problema
Quero ajustar os parâmetros do modelo de uma população simples de mistura 2-Gaussiana. Dado todo o hype em torno dos métodos bayesianos, quero entender se, para esse problema, a inferência bayesiana é uma ferramenta melhor que os métodos de ajuste tradicionais.
Até agora, o MCMC tem um desempenho muito ruim neste exemplo de brinquedo, mas talvez eu tenha esquecido alguma coisa. Então vamos ver o código.
As ferramentas
Vou usar python (2.7) + pilha scipy, lmfit 0.8 e PyMC 2.3.
Um caderno para reproduzir a análise pode ser encontrado aqui
Gere os dados
Primeiro vamos gerar os dados:
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])
O histograma de se samples
parece com isso:
um "pico amplo", os componentes são difíceis de detectar a olho nu.
Abordagem clássica: ajuste ao histograma
Vamos tentar a abordagem clássica primeiro. Usando o lmfit , é fácil definir um modelo de 2 picos:
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'
Finalmente, ajustamos o modelo ao algoritmo simplex:
fit_res = model.fit(data, x=x_data, method='nelder')
print fit_res.fit_report()
O resultado é a seguinte imagem (linhas tracejadas vermelhas são centros ajustados):
Mesmo que o problema seja meio difícil, com valores e restrições iniciais adequados, os modelos convergiram para uma estimativa bastante razoável.
Abordagem bayesiana: MCMC
Defino o modelo no PyMC de maneira hierárquica. centers
e sigmas
são as distribuições anteriores para os hiperparâmetros que representam os 2 centros e 2 sigmas dos 2 gaussianos. alpha
é a fração da primeira população e a distribuição anterior é aqui uma versão beta.
Uma variável categórica escolhe entre as duas populações. Entendo que essa variável precisa ter o mesmo tamanho dos dados ( samples
).
Finalmente mu
e tau
são deterministas variáveis que determinam os parâmetros da distribuição normal (que dependem da category
variável de modo que alternar aleatoriamente entre os dois valores para as duas populações).
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])
Então eu executo o MCMC com um número bastante longo de iterações (1e5, ~ 60s na minha máquina):
mcmc = pm.MCMC(model)
mcmc.sample(100000, 30000)
No entanto, os resultados são muito ímpares. Por exemplo, trace (a fração da primeira população) tende a 0 em vez de convergir para 0,4 e tem uma autocorrelação muito forte:
Também os centros dos gaussianos também não convergem. Por exemplo:
Como você vê na escolha anterior, tentei "ajudar" o algoritmo MCMC usando uma distribuição Beta para a fração populacional anterior . Também as distribuições anteriores para os centros e sigmas são bastante razoáveis (eu acho).
Então, o que está acontecendo aqui? Estou fazendo algo errado ou o MCMC não é adequado para esse problema?
Entendo que o método MCMC será mais lento, mas o ajuste trivial do histograma parece ter um desempenho imensamente melhor na resolução das populações.
fonte
proposal_distribution
eproposal_sd
e por que usarPrior
é melhor para as variáveis categóricas?Você não deve modelar com um Normal, dessa forma, permitindo valores negativos para a variação padrão. Em vez disso, use algo como:σ
atualizar:
Eu cheguei perto dos parâmetros iniciais dos dados alterando estas partes do seu modelo:
e invocando o mcmc com algum desbaste:
resultados:
Os posteriores não são muito agradáveis ... Na seção 11.6 do livro BUGS, eles discutem esse tipo de modelo e afirmam que existem problemas de convergência sem solução óbvia. Veja também aqui .
fonte
mcmc.sample(60000, 3000, 3)
Além disso, a não identificação é um grande problema para o uso do MCMC em modelos de mistura. Basicamente, se você alternar rótulos nos meios e nas atribuições do cluster, a probabilidade não muda e isso pode confundir o amostrador (entre cadeias ou dentro de cadeias). Uma coisa que você pode tentar diminuir é usar Potenciais no PyMC3. Uma boa implementação de um GMM com potenciais está aqui . O posterior para esses tipos de problemas também é geralmente altamente multimodal, o que também apresenta um grande problema. Você pode querer usar EM (ou inferência variacional).
fonte