Wie kann ich bei einer 10D-MCMC-Kette deren posterioren Modus (e) in R bestimmen?

10

Frage: Nehmen wir an, ich bin bereit, Ihnen mit einer 10-dimensionalen MCMC-Kette eine Matrix der Zeichnungen zu übergeben: 100.000 Iterationen (Zeilen) mit 10 Parametern (Spalten). Wie kann ich die posterioren Modi am besten identifizieren? Ich beschäftige mich besonders mit mehreren Modi.

Hintergrund:Ich betrachte mich als einen rechnerisch versierten Statistiker, aber als mir ein Kollege diese Frage stellte, schämte ich mich, dass ich keine vernünftige Antwort finden konnte. Das Hauptanliegen ist, dass mehrere Modi angezeigt werden können, jedoch nur, wenn mindestens acht der zehn Dimensionen berücksichtigt werden. Mein erster Gedanke wäre, eine Schätzung der Kerneldichte zu verwenden, aber eine Suche in R ergab nichts, was für Probleme mit mehr als drei Dimensionen vielversprechend wäre. Der Kollege hat eine Ad-hoc-Binning-Strategie in zehn Dimensionen vorgeschlagen und nach einem Maximum gesucht. Ich befürchte jedoch, dass die Bandbreite entweder zu erheblichen Sparsity-Problemen oder zu einer mangelnden Auflösung bei der Erkennung mehrerer Modi führen kann. Trotzdem würde ich gerne Vorschläge für automatisierte Bandbreitenvorschläge, Links zu einem 10-Kernel-Dichteschätzer oder alles andere, was Sie wissen, annehmen.

Sorgen:

  1. Wir glauben, dass die Verteilung ziemlich verzerrt sein kann; Daher möchten wir die posterioren Modi und nicht die posterioren Mittel identifizieren.

  2. Wir sind besorgt, dass es mehrere hintere Modi geben könnte.

  3. Wenn möglich, würden wir einen R-basierten Vorschlag bevorzugen. Aber jeder Algorithmus wird funktionieren, solange es nicht unglaublich schwierig ist, ihn zu implementieren. Ich denke, ich würde es vorziehen, keinen Nd-Kernel-Dichteschätzer mit automatisierter Bandbreitenauswahl von Grund auf neu zu implementieren.

M. Tibbits
quelle
Bitte lesen
Pavel Ruzankin

Antworten:

9

Haben Sie überlegt, einen Ansatz für den nächsten Nachbarn zu verwenden?

Erstellen Sie beispielsweise eine Liste der knächsten Nachbarn für jeden der 100'000 Punkte und betrachten Sie dann den Datenpunkt mit der kleinsten Entfernung des kthNachbarn als Modus. Mit anderen Worten: Finden Sie den Punkt mit der 'kleinsten Blase', die kandere Punkte um diesen Punkt enthält.

Ich bin mir nicht sicher, wie robust dies ist und die Wahl für kbeeinflusst offensichtlich die Ergebnisse.

Andre Holzner
quelle
Manchmal möchte ich mich nur auf den Kopf schlagen. Hervorragender Vorschlag.
M. Tibbits
1
Ich habe auch nur daran gedacht, die kmeansFunktion in R zu verwenden. Ich sollte zwischen Mitternacht und 4 Uhr morgens wirklich keine Fragen stellen.
M. Tibbits
4

Dies ist nur eine teilweise Antwort.

Ich habe kürzlich figtree für mehrdimensionale Kernel-Dichteschätzungen verwendet. Es ist ein C-Paket und ich habe es ziemlich einfach zum Laufen gebracht. Ich habe es jedoch nur verwendet, um die Dichte an bestimmten Punkten zu schätzen, nicht um zusammenfassende Statistiken zu berechnen.

csgillespie
quelle
3

Wenn Sie die Protokollwahrscheinlichkeiten beibehalten, können Sie einfach die mit dem höchsten Wert auswählen. Wenn Sie sich hauptsächlich für den Modus interessieren, reicht es aus, nur eine Optimierung durchzuführen, um den Punkt mit der höchsten Protokollwahrscheinlichkeit zu finden.

John Salvatier
quelle
Dies ist die relevanteste Antwort, zumindest der erste Teil! In vielen MCMC-Simulationen werden die (log-) Wahrscheinlichkeiten für alle Vorschläge berechnet und können somit gespeichert werden. Oder der bisher höchste Wert und sein Argument können gespeichert werden. Vorausgesetzt, der MCMC-Algorithmus hat über die Anzahl der von Ihnen ausgeführten Simulationen konvergiert, ist dies ein gültiger Ansatz.
Xi'an
2

Haben Sie über "PRIM / Bump Hunting" nachgedacht? (Siehe z. B. Abschnitt 9.3 von 'Die Elemente des statistischen Lernens' von Tibshirani et al. oder fragen Sie Ihre bevorzugte Suchmaschine). Ich bin mir nicht sicher, ob das in R implementiert ist.

[Soweit ich verstanden habe, versuchen Sie, den Modus der Wahrscheinlichkeitsdichte zu finden, aus der Ihre 100'000 Zeilen gezogen werden. Ihr Problem würde also teilweise gelöst, indem Sie eine geeignete density estimationMethode finden.

Andre Holzner
quelle
Ja, es gibt ein Prim- Paket mit einer R-Vignette: Verwenden von Prim für die Bump-Jagd . Mir ist jedoch nicht klar, wie es in diesem Fall funktionieren wird.
Chl