Würde eine "Wichtigkeits-Gibbs" -Sampling-Methode funktionieren?

8

Ich vermute, dass dies eine ziemlich ungewöhnliche und explorative Frage ist. Bitte nehmen Sie Kontakt mit mir auf.

Ich frage mich, ob man die Idee der Wichtigkeitsstichprobe auf Gibbs-Stichproben anwenden könnte. Ich meine Folgendes: Bei der Gibbs-Abtastung ändern wir jeweils den Wert einer Variablen (oder eines Variablenblocks), wobei wir aus der bedingten Wahrscheinlichkeit für die verbleibenden Variablen abtasten.

Es ist jedoch möglicherweise nicht möglich oder einfach, aus der genauen bedingten Wahrscheinlichkeit eine Stichprobe zu erstellen. Stattdessen probieren wir aus einer Angebotsverteilung und verwenden beispielsweise Metropolis-Hastings (MH).q

So weit, ist es gut. Aber hier ist ein abweichender Weg: Was passiert, wenn wir anstelle von MH dieselbe Idee verwenden, die bei der Wichtigkeitsstichprobe verwendet wird, nämlich aus und ein Wichtigkeitsgewicht der aktuellen Stichprobe beizubehalten?qp/q

Im Detail: Nehmen wir an, wir haben die Variablen und eine faktorisierte Verteilung so dass . Wir behalten die Vorschlagswahrscheinlichkeit bei, mit der der aktuelle Wert jeder Variablen wird . Bei jedem Schritt ändern wir eine Teilmenge der Variablen und aktualisieren (nur die Faktoren von und , die betroffen sind). Wir nehmen die Stichproben und ihr Wichtigkeitsgewicht, um die Statistik zu berechnen, an der wir interessiert sind.x1,,xnϕ1,,ϕmpi=1mϕiqixip(x)/q(x)pq

Wäre dieser Algorithmus korrekt? Wenn nicht, klare Gründe, warum nicht? Intuitiv macht es für mich Sinn, da es das gleiche zu tun scheint, was das Sampling von Wichtigkeit tut, aber stattdessen mit abhängigen Samples.

Ich habe dies für ein Gaußsches Random-Walk-Modell implementiert und festgestellt, dass die Gewichte immer kleiner (aber nicht monoton) werden, sodass die anfänglichen Stichproben zu wichtig sind und die Statistik dominieren. Ich bin mir ziemlich sicher, dass die Implementierung nicht fehlerhaft ist, da ich bei jedem Schritt das aktualisierte Gewicht mit einer expliziten Brute-Force-Berechnung vergleiche. Es ist zu beachten, dass die Gewichte nicht auf unbestimmte Zeit auf Null abfallen, da sie wobei sowohl als auch Produkte einer endlichen Anzahl von Dichten sind und jede Probe aus einer Normalverteilung erhalten wird, die nur selten Null ist.p/qpq

Ich versuche zu verstehen, warum die Gewichte so sinken und ob dies eine Folge davon ist, dass diese Methode tatsächlich nicht korrekt ist.


Hier ist eine genauere Definition des Algorithmus, der auf einen Gaußschen Zufallslauf für die Variablen angewendet wird . Der Code folgt unten.X1,,Xn

Das Modell ist einfach , wobei auf .XiN(Xi1,σ2),i=1,,nX00

Das Gewicht der aktuellen Stichprobe ist , wobei die Gaußschen Dichten und die Verteilungen sind, aus denen die aktuellen Werte abgetastet wurden. Zunächst werden die Werte einfach vorwärts abgetastet, sodass und das Anfangsgewicht .ip(xi)iq(xi)pqq=p1

Dann wähle ich bei jedem Schritt zum Ändern. Ich probiere einen neuen Wert für aus , damit diese Dichte zur neuen verwendeten Angebotsverteilung für .j{1,,n}xjXjN(Xj1,σ2)Xj

Um das Gewicht zu aktualisieren, dividiere ich es durch die Dichten und des alten Wertes gemäß und und multipliziere mit den Dichten und des neuen Wertes gemäß und . Dies aktualisiert den Zähler des Gewichts.p(xj|xj1)p(xj+1|xj)xjxj1xj+1p(xj|xj1)p(xj+1|xj)xjxj1xj+1p

Um den Nenner zu aktualisieren , multipliziere ich das Gewicht mit dem alten Vorschlag (und entferne ihn so vom Nenner) und dividiere ihn durch .qq(xj)q(xj)

(Da ich aus der auf zentrierten Normalen , ist immer gleich sodass sie sich aufheben und die Implementierung dies tut nicht wirklich benutzen).xjxj1q(xj)p(xj|xj1)

Wie ich bereits erwähnt habe, vergleiche ich im Code diese inkrementelle Gewichtsberechnung mit der tatsächlichen expliziten Berechnung, nur um sicherzugehen.


Hier ist der Code als Referenz.

println("Original sample: " + currentSample);
int flippedVariablesIndex = 1 + getRandom().nextInt(getVariables().size() - 1);
println("Flipping: " + flippedVariablesIndex);
double oldValue = getValue(currentSample, flippedVariablesIndex);
NormalDistribution normalFromBack = getNormalDistribution(getValue(currentSample, flippedVariablesIndex - 1));
double previousP = normalFromBack.density(oldValue);
double newValue = normalFromBack.sample();
currentSample.set(getVariable(flippedVariablesIndex), newValue);
double previousQ = fromVariableToQ.get(getVariable(flippedVariablesIndex));
fromVariableToQ.put(getVariable(flippedVariablesIndex), normalFromBack.density(newValue));
if (flippedVariablesIndex < length - 1) {
    NormalDistribution normal = getNormalDistribution(getValue(currentSample, flippedVariablesIndex + 1));
    double oldForwardPotential = normal.density(oldValue);
    double newForwardPotential = normal.density(newValue);
    // println("Removing old forward potential " + oldForwardPotential);
    currentSample.removePotential(new DoublePotential(oldForwardPotential));
    // println("Multiplying new forward potential " + newForwardPotential);
    currentSample.updatePotential(new DoublePotential(newForwardPotential));
}

// println("Removing old backward potential " + previousP);
currentSample.removePotential(new DoublePotential(previousP));
// println("Multiplying (removing from divisor) old q " + previousQ);
currentSample.updatePotential(new DoublePotential(previousQ));

println("Final sample: " + currentSample);
println();

// check by comparison to brute force calculation of weight:
double productOfPs = 1.0;
for (int i = 1; i != length; i++) {
    productOfPs *= getNormalDistribution(getValue(currentSample, i - 1)).density(getValue(currentSample, i));
}
double productOfQs = Util.fold(fromVariableToQ.values(), (p1, p2) -> p1*p2, 1.0);
double weight = productOfPs/productOfQs;
if (Math.abs(weight - currentSample.getPotential().doubleValue()) > 0.0000001) {
    println("Error in weight calculation");
    System.exit(0);
}
user118967
quelle
Die Wichtigkeitsabtastung liefert keine Stichproben aus der Zielverteilung (in diesem Fall die vollständigen Bedingungen von ). Die Markov-Kerneldynamik, die zu MCMC-Konvergenz führt, ist also nicht gültig. Ohne auf Ihren Code zu schauen, kann ich nicht sehen, warum die Gewichte auf 0 gehen.ϕi
Greenparker
Vielen Dank. Ich denke, ich muss mich mit den Theoremen der MCMC-Konvergenz befassen. Ich habe den Code für alle Fälle eingefügt, es ist ziemlich einfach. Vielen Dank.
user118967
1
Können Sie erklären, wie Sie den Algorithmus implementieren, anstatt den Rohcode (oder zusätzlich) einzuschließen? Was ist die Zielverteilung, was sind die vollständigen Bedingungen, was ist die Angebotsverteilung, wie kombinieren Sie die Gewichte usw. usw.
Greenparker
Vielen Dank. Ich habe es getan, bitte lassen Sie mich wissen, wenn dies irgendwo verwirrend ist.
user118967
@ Xi'an: Hier wird die Wichtigkeitsabtastung auf den Flip einer einzelnen Variablen angewendet. Anstatt den Vorschlag zu akzeptieren oder nicht wie in Metropolis Hastings, akzeptieren wir ihn immer, behalten aber ein Maß für die Wichtigkeit dieses Flip bei, indem wir die Wahrscheinlichkeit p durch den Vorschlag q für die zu spiegelnde Variable dividieren.
user118967

Antworten:

4

Dies ist eine interessante Idee, aber ich sehe einige Schwierigkeiten damit:

  1. Im Gegensatz zur Standard-Wichtigkeitsstichprobe oder sogar zur Metropolisierten Wichtigkeitsstichprobe wirkt der Vorschlag nicht im selben Raum wie die Zielverteilung, sondern in einem Raum mit kleinerer Dimension, sodass die Validierung unklar ist [und möglicherweise die Gewichte über Iterationen hinweg beibehalten muss, was zu einer Entartung führt].
  2. Die fehlenden Normalisierungskonstanten in den vollständigen Bedingungen ändern sich bei jeder Iteration, werden jedoch nicht berücksichtigt [siehe unten].
  3. Die Gewichte sind nicht begrenzt, da es entlang der Iterationen schließlich zu Simulationen mit einem sehr großen Gewicht kommen wird, es sei denn, man verfolgt das letzte Auftreten einer Aktualisierung für denselben Index , was mit der Markovschen Validierung des Gibbs-Samplers in Konflikt geraten kann . Das Ausführen eines bescheidenen Experiments mit und Iterationen zeigt einen Bereich von Gewichten von bis .jn=2T=1037.656397e-073.699364e+04

Um mehr Details zu erfahren, betrachten Sie ein zweidimensionales Ziel , einschließlich der richtigen Normalisierungskonstante, und implementieren Sie den Gibbs-Sampler mit den Vorschlägen und . Richtige Wichtigkeitsgewichte [im Sinne der Erzeugung der richtigen Erwartung, dh eines unverzerrten Schätzers für eine beliebige Funktion von ] für aufeinanderfolgende Simulationen sind entweder wobei und die von . Oder äquivalent p(,)qX(|y)qY(|x)(X,Y)

p(xt,yt1)qX(xt|yt1)mY(yt1)orp(xt1,yt)qY(yt|xt1)mX(xt1)
mX()mY()p(,)
pX(xt|yt1)qX(xt|yt1)orpY(yt|xt1)qY(yt|xt1)
In beiden Fällen erfordert dies die [unlösbaren] Randdichten von und unter dem Ziel .XYp(,)

Es lohnt sich zu vergleichen, was hier passiert, mit dem parallelen gewichteten Metropolis-Algorithmus . (Siehe zum Beispiel Schuster und Klebanov, 2018. ) Wenn das Ziel wieder und der Vorschlag , ist das Wichtigkeitsgewicht ist korrekt [zur Erstellung einer unvoreingenommenen Schätzung] und aktualisiert das frühere Gewicht nicht, sondern beginnt bei jeder Iteration von vorne.p(,)q(,|x,y)

p(x,y)q(x,y|x,y)

(C.) Eine Korrektur des ursprünglichen Gibbs-Vorschlags besteht darin, einen neuen Wert für den gesamten Vektor, z. B. , aus dem Gibbs-Vorschlag , weil dann das Wichtigkeitsgewicht korrekt ist [eine mögliche Normalisierung fehlt Konstante, die jetzt wirklich konstant ist und nicht aus früheren Gibbs-Iterationen stammt] .(x,y)qX(xt|yt1)qY(yt|xt)

p(xt,yt)qX(xt|yt1)qY(yt|xt)

Ein letzter Hinweis: Für das im Code berücksichtigte Random-Walk-Ziel ist eine direkte Simulation durch Kaskadierung möglich: Simulieren Sie , dann mit , & tc.X1X2X1

Xi'an
quelle