Le problème

Je veux ajuster les paramètres du modèle d'une population de mélange 2-gaussien simple. Étant donné tout le battage médiatique autour des méthodes bayésiennes, je veux comprendre si pour ce problème l'inférence bayésienne est un meilleur outil que les méthodes d'ajustement traditionnelles.

Jusqu'à présent, MCMC fonctionne très mal dans cet exemple de jouet, mais peut-être que j'ai juste oublié quelque chose. Voyons donc le code.

Les outils

J'utiliserai python (2.7) + pile scipy, lmfit 0.8 et PyMC 2.3.

Un cahier pour reproduire l'analyse peut être trouvé ici

Générez les données

Générons d'abord les données:

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

L'histogramme de samplesressemble à ceci:

un "pic large", les composants sont difficiles à repérer à l'œil nu.

Approche classique: ajuster l'histogramme

Essayons d'abord l'approche classique. En utilisant lmfit, il est facile de définir un modèle à 2 pics:

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'

Enfin, nous ajustons le modèle avec l'algorithme simplex:

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

Le résultat est l'image suivante (les lignes pointillées rouges sont des centres ajustés):

Même si le problème est assez difficile, avec des valeurs initiales et des contraintes appropriées, les modèles ont convergé vers une estimation assez raisonnable.

Approche bayésienne: MCMC

Je définis le modèle dans PyMC de manière hiérarchique. centerset sigmassont la distribution a priori des hyperparamètres représentant les 2 centres et 2 sigmas des 2 gaussiens. alphaest la fraction de la première population et la distribution antérieure est ici une Bêta.

Une variable catégorielle choisit entre les deux populations. Je crois comprendre que cette variable doit avoir la même taille que les données ( samples).

Enfin mu, ce tausont des variables déterministes qui déterminent les paramètres de la distribution normale (elles dépendent de la categoryvariable et basculent donc aléatoirement entre les deux valeurs pour les deux populations).

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

Ensuite, je lance le MCMC avec un nombre d'itérations assez long (1e5, ~ 60s sur ma machine):

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

$\alpha$

De plus, les centres des gaussiens ne convergent pas aussi bien. Par exemple:

$\alpha$

Alors qu'est-ce qui se passe ici? Suis-je en train de faire quelque chose de mal ou MCMC ne convient pas à ce problème?

Je comprends que la méthode MCMC sera plus lente, mais l'ajustement de l'histogramme trivial semble être extrêmement efficace pour résoudre les populations.

Le problème est dû à la façon dont PyMC tire des échantillons pour ce modèle. Comme expliqué dans la section 5.8.1 de la documentation PyMC, tous les éléments d'une variable de tableau sont mis à jour ensemble. Pour les petits tableaux comme centercelui-ci, ce n'est pas un problème, mais pour un grand tableau comme categorycelui-ci, cela conduit à un faible taux d'acceptation. Vous pouvez voir le taux d'acceptation via

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

La solution suggérée dans la documentation consiste à utiliser un tableau de variables à valeurs scalaires. De plus, vous devez configurer certaines des distributions de proposition car les choix par défaut sont mauvais. Voici le code qui me convient:

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

Les options proposal_sdet proposal_distributionsont expliquées dans la section 5.7.1 . Pour les centres, j'ai défini la proposition pour correspondre à peu près à l'écart-type de la partie postérieure, qui est beaucoup plus petit que la valeur par défaut en raison de la quantité de données. PyMC essaie de régler la largeur de la proposition, mais cela ne fonctionne que si votre taux d'acceptation est suffisamment élevé pour commencer. Pour category, la valeur par défaut proposal_distribution = 'Poisson'qui donne de mauvais résultats (je ne sais pas pourquoi, mais cela ne ressemble certainement pas à une proposition sensée pour une variable binaire).

$\sigma$

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

mettre à jour:

Je me suis approché des paramètres initiaux des données en modifiant ces parties de votre modèle:

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

et en invoquant le mcmc avec un éclaircissement:

mcmc.sample(200000, 3000, 10)

résultats:

Les postérieurs ne sont pas très gentils ... Dans la section 11.6 du livre BUGS, ils discutent de ce type de modèle et déclarent qu'il existe des problèmes de convergence sans solution évidente. Vérifiez également ici .

En outre, la non-identifiabilité est un gros problème pour l'utilisation de MCMC pour les modèles de mélange. Fondamentalement, si vous changez d'étiquettes sur vos moyennes de cluster et les affectations de cluster, la probabilité ne change pas, ce qui peut confondre l'échantillonneur (entre les chaînes ou au sein des chaînes). Une chose que vous pourriez essayer d'atténuer est d'utiliser les potentiels dans PyMC3. Une bonne implémentation d'un GMM avec des potentiels est ici . La partie postérieure de ces types de problèmes est également généralement très multimodale, ce qui présente également un gros problème. Vous pouvez plutôt utiliser EM (ou l'inférence variationnelle).