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 samples
sieht folgendermaßen aus:
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):
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. centers
und sigmas
sind die Prioritätenverteilung für die Hyperparameter, die die 2 Zentren und 2 Sigmen der 2 Gaußschen darstellen. alpha
ist 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 mu
und tau
deterministisch sind Variablen, die die Parameter der Normalverteilung zu bestimmen (sie sind abhängig von den category
variablen , 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:
Auch die Zentren der Gaußschen konvergieren nicht so gut. Zum Beispiel:
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.
quelle
proposal_distribution
undproposal_sd
warum verwendenPrior
ist besser für die kategorischen Variablen?aktualisieren:
Ich habe mich den Anfangsparametern der Daten genähert, indem ich diese Teile Ihres Modells geändert habe:
und durch Aufrufen des mcmc mit etwas Ausdünnung:
Ergebnisse:
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 .
quelle
mcmc.sample(60000, 3000, 3)
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.
quelle