2-Gaußsche Mischungsmodellinferenz mit MCMC und PyMC

10

Das Problem

Ich möchte die Modellparameter einer einfachen 2-Gaußschen Mischungspopulation anpassen. Angesichts des Hype um Bayes'sche Methoden möchte ich verstehen, ob die Bayes'sche Inferenz für dieses Problem ein besseres Werkzeug ist als herkömmliche Anpassungsmethoden.

Bisher ist MCMC in diesem Spielzeugbeispiel sehr schlecht, aber vielleicht habe ich gerade etwas übersehen. Schauen wir uns also den Code an.

Die Werkzeuge

Ich werde Python (2.7) + Scipy Stack, lmfit 0.8 und PyMC 2.3 verwenden.

Ein Notizbuch zur Reproduktion der Analyse finden Sie hier

Generieren Sie die Daten

Lassen Sie uns zuerst die Daten generieren:

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

Das Histogramm von samplessieht folgendermaßen aus:

Datenhistogramm

Bei einem "breiten Peak" sind die Komponenten mit dem Auge schwer zu erkennen.

Klassischer Ansatz: Passen Sie das Histogramm an

Versuchen wir zuerst den klassischen Ansatz. Mit lmfit ist es einfach, ein 2-Peaks-Modell zu definieren:

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'

Schließlich passen wir das Modell mit dem Simplex-Algorithmus an:

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

Das Ergebnis ist das folgende Bild (rote gestrichelte Linien sind angepasste Zentren):

NLS-Fit-Ergebnisse

Selbst wenn das Problem schwierig ist, konvergierten die Modelle mit den richtigen Anfangswerten und Einschränkungen zu einer vernünftigen Schätzung.

Bayesianischer Ansatz: MCMC

Ich definiere das Modell in PyMC hierarchisch. centersund sigmassind die Prioritätenverteilung für die Hyperparameter, die die 2 Zentren und 2 Sigmen der 2 Gaußschen darstellen. alphaist der Anteil der ersten Population und die vorherige Verteilung ist hier eine Beta.

Eine kategoriale Variable wählt zwischen den beiden Populationen. Nach meinem Verständnis muss diese Variable dieselbe Größe wie die Daten haben ( samples).

Schließlich muund taudeterministisch sind Variablen, die die Parameter der Normalverteilung zu bestimmen (sie sind abhängig von den categoryvariablen , so dass sie zufällig zwischen den beiden Werten für die zwei Populationen Schalter).

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

Dann führe ich die MCMC mit einer ziemlich langen Anzahl von Iterationen aus (1e5, ~ 60s auf meinem Computer):

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

Die Ergebnisse sind jedoch sehr merkwürdig. Zum Beispiel tendiert trace (der Anteil der ersten Population) gegen 0, um gegen 0,4 zu konvergieren, und hat eine sehr starke Autokorrelation:α

MCMC Alpha-Zusammenfassung

Auch die Zentren der Gaußschen konvergieren nicht so gut. Zum Beispiel:

Zusammenfassung der MCMC-Zentren_0

Wie Sie in der vorherigen Auswahl sehen, habe ich versucht, den MCMC-Algorithmus mithilfe einer Beta-Verteilung für die vorherige Populationsfraktion zu "unterstützen" . Auch die vorherigen Verteilungen für die Zentren und Sigmen sind ziemlich vernünftig (denke ich).α

Also, was ist hier los? Mache ich etwas falsch oder ist MCMC für dieses Problem nicht geeignet?

Ich verstehe, dass die MCMC-Methode langsamer sein wird, aber die triviale Histogrammanpassung scheint bei der Auflösung der Populationen immens besser zu funktionieren.

user2304916
quelle

Antworten:

6

Das Problem wird durch die Art und Weise verursacht, wie PyMC Proben für dieses Modell zeichnet. Wie in Abschnitt 5.8.1 der PyMC-Dokumentation erläutert, werden alle Elemente einer Array-Variablen zusammen aktualisiert. Für kleine Arrays wie centerdieses ist dies kein Problem, aber für ein großes Array wie categorydieses führt dies zu einer geringen Akzeptanzrate. Sie können die Akzeptanzrate über sehen

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

Die in der Dokumentation vorgeschlagene Lösung besteht darin, ein Array skalarwertiger Variablen zu verwenden. Darüber hinaus müssen Sie einige der Angebotsverteilungen konfigurieren, da die Standardauswahl schlecht ist. Hier ist der Code, der für mich funktioniert:

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

Die proposal_sdund proposal_distributionOptionen werden in Abschnitt 5.7.1 erläutert . Für die Zentren habe ich den Vorschlag so festgelegt, dass er ungefähr der Standardabweichung des Seitenzahns entspricht, die aufgrund der Datenmenge viel kleiner als die Standardabweichung ist. PyMC versucht zwar, die Breite des Vorschlags anzupassen, dies funktioniert jedoch nur, wenn Ihre Akzeptanzrate zunächst hoch genug ist. Denn categorydie Standardeinstellung, proposal_distribution = 'Poisson'die schlechte Ergebnisse liefert (ich weiß nicht, warum das so ist, aber es klingt sicherlich nicht nach einem vernünftigen Vorschlag für eine binäre Variable).

Tom Minka
quelle
Danke, das ist wirklich nützlich, obwohl es fast unerträglich langsam wird. Können Sie kurz über die Bedeutung von erklären proposal_distributionund proposal_sdwarum verwenden Priorist besser für die kategorischen Variablen?
user2304916
Danke, Tom. Ich stimme zu, dass Poisson hier eine seltsame Wahl ist. Ich habe ein Problem eröffnet: github.com/pymc-devs/pymc/issues/627
twiecki
2

σ

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

aktualisieren:

Ich habe mich den Anfangsparametern der Daten genähert, indem ich diese Teile Ihres Modells geändert habe:

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

und durch Aufrufen des mcmc mit etwas Ausdünnung:

mcmc.sample(200000, 3000, 10)

Ergebnisse:

Alpha

Zentren

Sigmas

Die Posterioren sind nicht sehr nett, du ... In Abschnitt 11.6 des BUGS-Buches diskutieren sie diese Art von Modell und geben an, dass es Konvergenzprobleme gibt, für die es keine offensichtliche Lösung gibt. Überprüfen Sie auch hier .

jpneto
quelle
Das ist ein guter Punkt, ich benutze jetzt ein Gamma. Trotzdem tendiert die Alpha-Kurve immer zu 0 (anstelle von 0,4). Ich frage mich, ob in meinem Beispiel ein dummer Fehler lauert.
user2304916
Ich habe Gamma (.1, .1) ausprobiert, aber Exp (.1) scheint besser zu funktionieren. Wenn die Autokorrelation hoch ist, können Sie auch eine gewisse Ausdünnung einschließen, z. B.mcmc.sample(60000, 3000, 3)
jpneto
2

Die Nichtidentifizierbarkeit ist auch ein großes Problem bei der Verwendung von MCMC für Mischungsmodelle. Grundsätzlich ändert sich die Wahrscheinlichkeit nicht, wenn Sie die Beschriftungen Ihrer Cluster-Mittelwerte und Cluster-Zuweisungen wechseln. Dies kann den Sampler (zwischen Ketten oder innerhalb von Ketten) verwirren. Eine Möglichkeit , dies zu mildern, ist die Verwendung von Potenzialen in PyMC3. Eine gute Implementierung eines GMM mit Potenzialen ist hier . Der hintere Teil für diese Art von Problemen ist im Allgemeinen auch sehr multimodal, was ebenfalls ein großes Problem darstellt. Möglicherweise möchten Sie stattdessen EM (oder Variationsinferenz) verwenden.

Benjamin Doughty
quelle