El problema

Quiero ajustar los parámetros del modelo de una población simple de mezcla 2-Gaussiana. Dada toda la exageración en torno a los métodos bayesianos, quiero entender si para este problema la inferencia bayesiana es una herramienta mejor que los métodos de ajuste tradicionales.

Hasta ahora, MCMC funciona muy mal en este ejemplo de juguete, pero tal vez simplemente pasé por alto algo. Entonces, veamos el código.

Las herramientas

Usaré python (2.7) + scipy stack, lmfit 0.8 y PyMC 2.3.

Puede encontrar un cuaderno para reproducir el análisis aquí.

Generar los datos

Primero vamos a generar los datos:

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])

El histograma de se samplesve así:

un "pico amplio", los componentes son difíciles de detectar a simple vista.

Enfoque clásico: ajustar el histograma

Probemos el enfoque clásico primero. Usando lmfit es fácil definir un 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 el modelo con el algoritmo simplex:

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

El resultado es la siguiente imagen (las líneas discontinuas rojas son centros ajustados):

Incluso si el problema es un poco difícil, con valores iniciales adecuados y restricciones, los modelos convergieron en una estimación bastante razonable.

Enfoque bayesiano: MCMC

Defino el modelo en PyMC de forma jerárquica. centersy sigmasson la distribución previa de los hiperparámetros que representan los 2 centros y 2 sigmas de los 2 gaussianos. alphaes la fracción de la primera población y la distribución anterior es aquí una Beta.

Una variable categórica elige entre las dos poblaciones. Entiendo que esta variable debe tener el mismo tamaño que los datos ( samples).

Finalmente muy tauson variables deterministas que determinan los parámetros de la distribución Normal (dependen de la categoryvariable por lo que cambian aleatoriamente entre los dos valores para las dos poblaciones).

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])

Luego ejecuto el MCMC con un número bastante largo de iteraciones (1e5, ~ 60s en mi máquina):

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

$\alpha$

Además, los centros de los gaussianos no convergen también. Por ejemplo:

$\alpha$

Entonces, ¿qué está pasando aquí? ¿Estoy haciendo algo mal o MCMC no es adecuado para este problema?

Entiendo que el método MCMC será más lento, pero el ajuste del histograma trivial parece funcionar inmensamente mejor en la resolución de las poblaciones.

El problema es causado por la forma en que PyMC dibuja muestras para este modelo. Como se explica en la sección 5.8.1 de la documentación de PyMC, todos los elementos de una variable de matriz se actualizan juntos. Para arreglos pequeños como centereste no es un problema, pero para arreglos grandes como categoryesto conduce a una baja tasa de aceptación. Puede ver la tasa de aceptación a través de

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

La solución sugerida en la documentación es utilizar una matriz de variables con valores escalares. Además, debe configurar algunas de las distribuciones de propuestas, ya que las opciones predeterminadas son malas. Aquí está el código que funciona para mí:

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

Las opciones proposal_sdy proposal_distributionse explican en la sección 5.7.1 . Para los centros, configuré la propuesta para que coincida aproximadamente con la desviación estándar de la parte posterior, que es mucho más pequeña que la predeterminada debido a la cantidad de datos. PyMC intenta ajustar el ancho de la propuesta, pero esto solo funciona si su tasa de aceptación es lo suficientemente alta para empezar. Para category, el valor predeterminado proposal_distribution = 'Poisson'que da malos resultados (no sé por qué es esto, pero ciertamente no parece una propuesta sensata para una variable binaria).

$\sigma$

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

actualizar:

Me acerqué a los parámetros iniciales de los datos cambiando estas partes de su modelo:

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

e invocando el mcmc con algo de adelgazamiento:

mcmc.sample(200000, 3000, 10)

resultados:

Los posteriores no son muy agradables ... En la sección 11.6 del Libro BUGS , discuten este tipo de modelo y afirman que hay problemas de convergencia sin una solución obvia. Consulta también aquí .

Además, la no identificabilidad es un gran problema para usar MCMC para modelos mixtos. Básicamente, si cambia las etiquetas en sus medios de clúster y las asignaciones de clúster, la probabilidad no cambia, y esto puede confundir al muestreador (entre cadenas o dentro de cadenas). Una cosa que podría intentar mitigar esto es usar los potenciales en PyMC3. Una buena implementación de un GMM con potenciales está aquí . El posterior para este tipo de problemas también es generalmente altamente multimodal, lo que también presenta un gran problema. Es posible que desee utilizar EM (o inferencia variacional).