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

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?$q$$p/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.$x_1,\dots,x_n$$\phi_1,\dots,\phi_m$$p \propto \prod_{i=1}^m \phi_i$$q_i$$x_i$$p(x)/q(x)$$p$$q$

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/q$$p$$q$

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.$X_1,\dots,X_n$

Das Modell ist einfach , wobei auf .$X_i \sim \mathcal N(X_{i-1}, \sigma^2), i = 1,\dots,n$$X_0$$0$

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 .$\frac{\prod_i p(x_i)}{\prod_i q(x_i)}$$p$$q$$q = p$$1$

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 \in \{1,\dots,n\}$$x'_j$$X_j$$\mathcal N(X_{j-1},\sigma^2)$$X_j$

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(x_j | x_{j-1})$$p(x_{j+1} | x_j)$$x_j$$x_{j-1}$$x_{j+1}$$p(x'_j | x_{j-1})$$p(x_{j+1} | x'_j)$$x'_j$$x_{j-1}$$x_{j+1}$$p$

Um den Nenner zu aktualisieren , multipliziere ich das Gewicht mit dem alten Vorschlag (und entferne ihn so vom Nenner) und dividiere ihn durch .$q$$q(x_j)$$q(x'_j)$

(Da ich aus der auf zentrierten Normalen , ist immer gleich sodass sie sich aufheben und die Implementierung dies tut nicht wirklich benutzen).$x'_j$$x_{j-1}$$q(x'_j)$$p(x'_j | x_{j-1})$

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);
}

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 .$j$$n=2$$T=10^3$7.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(\cdot,\cdot)$$q_X(\cdot|y)$$q_Y(\cdot|x)$$(X,Y)$

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(\cdot,\cdot)$$q(\cdot,\cdot|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)$$q_X(x_t|y_{t-1})q_Y(y_t|x_{t})$

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.$X_1$$X_2$$X_1$