Verständnis von Metropolis-Hastings mit asymmetrischer Angebotsverteilung

Ich habe versucht, den Metropolis-Hastings-Algorithmus zu verstehen, um einen Code zum Schätzen der Parameter eines Modells zu schreiben (dh $f(x)=a*x$ ). Laut Bibliographie hat der Metropolis-Hastings-Algorithmus die folgenden Schritte:

• Generiere $Y_t \sim q(y|x^t)$
• $X^{t+1}=\begin{cases} Y^t, & \text{with probability} \quad \rho(x^t,Y_t), \\ x^t, & \text{with probability} \quad 1-\rho(x^t,Y_t), \end{cases}$

wobei $\rho(x,y)=\min \left( \frac{f(y)}{f(x)}*\frac{q(x|y)}{q(y|x)},1 \right)$

Wie möchte ich ein paar Fragen stellen:

• Die Bibliographie besagt, dass wenn eine symmetrische Verteilung ist, das Verhältnis q ( x | y ) / q ( y | x ) 1 wird und der Algorithmus Metropolis heißt. Ist das korrekt? Der einzige Unterschied zwischen Metropolis und Metropolis-Hastings besteht darin, dass die erste eine symmetrische Verteilung verwendet. Was ist mit "Random Walk" Metropolis (-Hastings)? Wie unterscheidet es sich von den anderen beiden?$q$$q(x|y)/q(y|x)$
• Der größte Teil des Beispielcodes, den ich online finde, verwendet eine Gaußsche Angebotsverteilung und daher ist ρ ( x , y ) = min ( f ( y ) / f ( x ) , 1 ) wobei f ( y ) / f ( x $q$$\rho(x,y)=\min( f(y)/f(x),1)$$f(y)/f(x)$ist das Wahrscheinlichkeitsverhältnis. Was ist, wenn die Angebotsverteilung eine Poisson-Verteilung ist? Ich denke, ich verstehe rational, warum dieses Verhältnis bei Verwendung einer asymmetrischen Verteilung nicht zu 1 wird, aber ich bin nicht ganz sicher, ob ich es mathematisch verstehe oder wie ich es mit Code implementieren soll. Könnte mir jemand einen einfachen Code (C, Python, R, Pseudo-Code oder was auch immer Sie bevorzugen) für den Metropolis-Hastings-Algorithmus mit einer nicht symmetrischen Angebotsverteilung geben?

Die Bibliographie besagt, dass wenn q eine symmetrische Verteilung ist, das Verhältnis q (x | y) / q (y | x) 1 wird und der Algorithmus Metropolis heißt. Ist das korrekt?

Ja das ist korrekt. Der Metropolis-Algorithmus ist ein Sonderfall des MH-Algorithmus.

Was ist mit "Random Walk" Metropolis (-Hastings)? Wie unterscheidet es sich von den anderen beiden?

Bei einem zufälligen Durchgang wird die Angebotsverteilung nach jedem Schritt auf den Wert neu zentriert, der zuletzt von der Kette generiert wurde. Im Allgemeinen ist die Angebotsverteilung bei einem zufälligen Spaziergang Gauß'sch. In diesem Fall erfüllt dieser zufällige Spaziergang die Symmetrieanforderung und der Algorithmus ist Metropolis. Ich nehme an, Sie könnten einen "Pseudo" -Zufallsrundgang mit einer asymmetrischen Verteilung durchführen, der dazu führen würde, dass die Vorschläge zu stark in die entgegengesetzte Richtung des Versatzes driften (eine linksversetzte Verteilung würde die Vorschläge nach rechts begünstigen). Ich bin nicht sicher, warum Sie dies tun würden, aber Sie könnten und es wäre ein Metropolis-Hastings-Algorithmus (dh erfordern Sie den zusätzlichen Verhältnis-Ausdruck).

Wie unterscheidet es sich von den anderen beiden?

In einem nicht zufälligen Walk-Algorithmus sind die Angebotsverteilungen festgelegt. In der Zufallsvariante ändert sich das Zentrum der Angebotsverteilung bei jeder Iteration.

Was ist, wenn die Angebotsverteilung eine Poisson-Verteilung ist?

Dann müssen Sie MH anstelle von Metropolis verwenden. Vermutlich wäre dies ein Beispiel für eine diskrete Verteilung, ansonsten würden Sie keine diskrete Funktion verwenden, um Ihre Vorschläge zu generieren.

In jedem Fall möchten Sie, wenn die Stichprobenverteilung abgeschnitten ist oder Sie die Abweichung bereits kennen, wahrscheinlich eine asymmetrische Stichprobenverteilung verwenden und müssen daher Metropolen-Hastings verwenden.

Könnte mir jemand ein einfaches Code-Beispiel geben (C, Python, R, Pseudo-Code oder was auch immer Sie bevorzugen)?

Hier ist die Metropole:

Metropolis <- function(F_sample # distribution we want to sample
                      , F_prop  # proposal distribution 
                      , I=1e5   # iterations
               ){
  y = rep(NA,T)
  y[1] = 0    # starting location for random walk
  accepted = c(1)

  for(t in 2:I)    {
    #y.prop <- rnorm(1, y[t-1], sqrt(sigma2) ) # random walk proposal
    y.prop <- F_prop(y[t-1]) # implementation assumes a random walk. 
                             # discard this input for a fixed proposal distribution

    # We work with the log-likelihoods for numeric stability.
    logR = sum(log(F_sample(y.prop))) -
           sum(log(F_sample(y[t-1])))    

    R = exp(logR)

    u <- runif(1)        ## uniform variable to determine acceptance
    if(u < R){           ## accept the new value
      y[t] = y.prop
      accepted = c(accepted,1)
    }    
    else{
      y[t] = y[t-1]      ## reject the new value
      accepted = c(accepted,0)
    }    
  }
  return(list(y, accepted))
}

Lassen Sie uns versuchen, mit diesem Beispiel eine bimodale Verteilung zu erstellen. Lassen Sie uns zuerst sehen, was passiert, wenn wir einen zufälligen Spaziergang für unseren Propsal verwenden:

set.seed(100)

test = function(x){dnorm(x,-5,1)+dnorm(x,7,3)}

# random walk
response1 <- Metropolis(F_sample = test
                       ,F_prop = function(x){rnorm(1, x, sqrt(0.5) )}
                      ,I=1e5
                       )
y_trace1 = response1[[1]]; accpt_1 = response1[[2]]
mean(accpt_1) # acceptance rate without considering burn-in
# 0.85585   not bad

# looks about how we'd expect
plot(density(y_trace1))
abline(v=-5);abline(v=7) # Highlight the approximate modes of the true distribution

Versuchen wir nun, Stichproben mit einer festen Angebotsverteilung zu erstellen und sehen, was passiert:

response2 <- Metropolis(F_sample = test
                            ,F_prop = function(x){rnorm(1, -5, sqrt(0.5) )}
                            ,I=1e5
                       )

y_trace2 = response2[[1]]; accpt_2 = response2[[2]]
mean(accpt_2) # .871, not bad

Das sieht auf den ersten Blick in Ordnung aus, aber wenn wir uns die posteriore Dichte ansehen ...

plot(density(y_trace2))

wir werden sehen, dass es bei einem lokalen Maximum vollständig stecken bleibt. Dies ist nicht ganz überraschend, da wir unsere Angebotsverteilung dort tatsächlich zentriert haben. Dasselbe passiert, wenn wir dies auf den anderen Modus zentrieren:

response2b <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, 7, sqrt(10) )}
                        ,I=1e5
)

plot(density(response2b[[1]]))

Wir können versuchen, unseren Vorschlag zwischen den beiden Modi fallen zu lassen, aber wir müssen die Varianz wirklich hoch einstellen, um eine Chance zu haben, einen von beiden zu erkunden

response3 <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, -2, sqrt(10) )}
                        ,I=1e5
)
y_trace3 = response3[[1]]; accpt_3 = response3[[2]]
mean(accpt_3) # .3958! 

Beachten Sie, dass die Auswahl des Zentrums unserer Angebotsverteilung einen erheblichen Einfluss auf die Akzeptanzrate unseres Samplers hat.

plot(density(y_trace3))

plot(y_trace3) # we really need to set the variance pretty high to catch 
               # the mode at +7. We're still just barely exploring it

Wir bleiben immer noch im engeren der beiden Modi stecken. Lassen Sie uns versuchen, dies direkt zwischen den beiden Modi abzulegen.

response4 <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, 1, sqrt(10) )}
                        ,I=1e5
)
y_trace4 = response4[[1]]; accpt_4 = response4[[2]]

plot(density(y_trace1))
lines(density(y_trace4), col='red')

Schließlich nähern wir uns dem, wonach wir gesucht haben. Wenn wir den Sampler lange genug laufen lassen, können wir theoretisch eine repräsentative Stichprobe aus jeder dieser Angebotsverteilungen erhalten, aber der Zufallsrundgang ergab sehr schnell eine brauchbare Stichprobe, und wir mussten unser Wissen darüber, wie der Posterior angenommen wurde, nutzen um die festen Stichprobenverteilungen auf ein brauchbares Ergebnis abzustimmen (was wir allerdings noch nicht ganz wissen)y_trace4 ).

Ich versuche es später mit einem Beispiel für Metropolen-Hastings. Sie sollten ziemlich leicht verstehen, wie Sie den obigen Code ändern können, um einen Metropolis-Hastings-Algorithmus zu erstellen (Hinweis: Sie müssen nur das zusätzliche Verhältnis zur logRBerechnung hinzufügen ).