Verwechselt mit MCMC Metropolis-Hastings-Variationen: Random-Walk, Non-Random-Walk, Independent, Metropolis

In den letzten Wochen habe ich versucht, MCMC und die Metropolis-Hastings-Algorithmen zu verstehen. Jedes Mal, wenn ich denke, dass ich es verstehe, stelle ich fest, dass ich falsch liege. Die meisten Codebeispiele, die ich online finde, implementieren etwas, das nicht mit der Beschreibung übereinstimmt. Dh: Sie sagen, sie implementieren Metropolis-Hastings, aber sie implementieren tatsächlich Metropolis mit wahlfreiem Zugriff. Andere (fast immer) überspringen stillschweigend die Implementierung des Hastings-Korrekturverhältnisses, weil sie eine symmetrische Angebotsverteilung verwenden. Eigentlich habe ich noch kein einfaches Beispiel gefunden, das das Verhältnis berechnet. Das verwirrt mich noch mehr. Kann mir jemand Codebeispiele (in einer beliebigen Sprache) für Folgendes geben:

• Vanille-Nicht-Zufalls-Metropolis-Hastings-Algorithmus mit Hastings-Korrekturverhältnisberechnung (auch wenn dies bei Verwendung einer symmetrischen Angebotsverteilung 1 ergibt).
• Vanilla Random Walk Metropolis-Hastings-Algorithmus.
• Vanilla Independent Metropolis-Hastings-Algorithmus.

Die Metropolis-Algorithmen müssen nicht angegeben werden, denn wenn ich mich nicht irre, besteht der einzige Unterschied zwischen Metropolis und Metropolis-Hastings darin, dass die ersten immer aus einer symmetrischen Verteilung abtasten und daher nicht das Hastings-Korrekturverhältnis aufweisen. Es ist keine detaillierte Erläuterung der Algorithmen erforderlich. Ich verstehe die Grundlagen, aber ich bin ein bisschen verwirrt mit all den verschiedenen Namen für die verschiedenen Variationen des Metropolis-Hastings-Algorithmus, aber auch damit, wie Sie das Hastings-Korrekturverhältnis auf dem Vanilla-Non-Random-Walk-MH praktisch implementieren. Bitte kopieren Sie keine Links zum Einfügen, die meine Fragen teilweise beantworten, da ich sie höchstwahrscheinlich bereits gesehen habe. Diese Links führten mich zu dieser Verwirrung. Vielen Dank.

Hier gehts - drei Beispiele. Ich habe den Code viel weniger effizient gemacht als in einer echten Anwendung, um die Logik klarer zu machen (hoffe ich).

# We'll assume estimation of a Poisson mean as a function of x
x <- runif(100)
y <- rpois(100,5*x)  # beta = 5 where mean(y[i]) = beta*x[i]

# Prior distribution on log(beta): t(5) with mean 2 
# (Very spread out on original scale; median = 7.4, roughly)
log_prior <- function(log_beta) dt(log_beta-2, 5, log=TRUE)

# Log likelihood
log_lik <- function(log_beta, y, x) sum(dpois(y, exp(log_beta)*x, log=TRUE))

# Random Walk Metropolis-Hastings 
# Proposal is centered at the current value of the parameter

rw_proposal <- function(current) rnorm(1, current, 0.25)
rw_p_proposal_given_current <- function(proposal, current) dnorm(proposal, current, 0.25, log=TRUE)
rw_p_current_given_proposal <- function(current, proposal) dnorm(current, proposal, 0.25, log=TRUE)

rw_alpha <- function(proposal, current) {
   # Due to the structure of the rw proposal distribution, the rw_p_proposal_given_current and
   # rw_p_current_given_proposal terms cancel out, so we don't need to include them - although
   # logically they are still there:  p(prop|curr) = p(curr|prop) for all curr, prop
   exp(log_lik(proposal, y, x) + log_prior(proposal) - log_lik(current, y, x) - log_prior(current))
}

# Independent Metropolis-Hastings
# Note: the proposal is independent of the current value (hence the name), but I maintain the
# parameterization of the functions anyway.  The proposal is not ignorable any more
# when calculation the acceptance probability, as p(curr|prop) != p(prop|curr) in general.

ind_proposal <- function(current) rnorm(1, 2, 1) 
ind_p_proposal_given_current <- function(proposal, current) dnorm(proposal, 2, 1, log=TRUE)
ind_p_current_given_proposal <- function(current, proposal) dnorm(current, 2, 1, log=TRUE)

ind_alpha <- function(proposal, current) {
   exp(log_lik(proposal, y, x)  + log_prior(proposal) + ind_p_current_given_proposal(current, proposal) 
       - log_lik(current, y, x) - log_prior(current) - ind_p_proposal_given_current(proposal, current))
}

# Vanilla Metropolis-Hastings - the independence sampler would do here, but I'll add something
# else for the proposal distribution; a Normal(current, 0.1+abs(current)/5) - symmetric but with a different
# scale depending upon location, so can't ignore the proposal distribution when calculating alpha as
# p(prop|curr) != p(curr|prop) in general

van_proposal <- function(current) rnorm(1, current, 0.1+abs(current)/5)
van_p_proposal_given_current <- function(proposal, current) dnorm(proposal, current, 0.1+abs(current)/5, log=TRUE)
van_p_current_given_proposal <- function(current, proposal) dnorm(current, proposal, 0.1+abs(proposal)/5, log=TRUE)

van_alpha <- function(proposal, current) {
   exp(log_lik(proposal, y, x)  + log_prior(proposal) + ind_p_current_given_proposal(current, proposal) 
       - log_lik(current, y, x) - log_prior(current) - ind_p_proposal_given_current(proposal, current))
}


# Generate the chain
values <- rep(0, 10000) 
u <- runif(length(values))
naccept <- 0
current <- 1  # Initial value
propfunc <- van_proposal  # Substitute ind_proposal or rw_proposal here
alphafunc <- van_alpha    # Substitute ind_alpha or rw_alpha here
for (i in 1:length(values)) {
   proposal <- propfunc(current)
   alpha <- alphafunc(proposal, current)
   if (u[i] < alpha) {
      values[i] <- exp(proposal)
      current <- proposal
      naccept <- naccept + 1
   } else {
      values[i] <- exp(current)
   }
}
naccept / length(values)
summary(values)

Für den Vanille-Sampler erhalten wir:

> naccept / length(values)
[1] 0.1737
> summary(values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.843   5.153   5.388   5.378   5.594   6.628 

Das ist eine geringe Akzeptanzwahrscheinlichkeit, aber dennoch ... würde es helfen, den Vorschlag zu optimieren oder einen anderen anzunehmen. Hier sind die Ergebnisse des zufälligen Wandervorschlags:

> naccept / length(values)
[1] 0.2902
> summary(values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.718   5.147   5.369   5.370   5.584   6.781 

Ähnliche Ergebnisse, wie man hoffen würde, und eine bessere Akzeptanzwahrscheinlichkeit (angestrebt werden ~ 50% mit einem Parameter.)

Und der Vollständigkeit halber der Independence Sampler:

> naccept / length(values)
[1] 0.0684
> summary(values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.990   5.162   5.391   5.380   5.577   8.802 

Da es sich nicht an die Form des Seitenzahns "anpasst", weist es tendenziell die geringste Akzeptanzwahrscheinlichkeit auf und ist für dieses Problem am schwierigsten zu optimieren.

Beachten Sie, dass wir im Allgemeinen Vorschläge mit dickeren Schwänzen bevorzugen, aber das ist ein ganz anderes Thema.

Sehen:

Konstruktionsbedingt hängt der Algorithmus nicht von der Normalisierungskonstante ab, da das Verhältnis der PDFs von Bedeutung ist. Die Variation des Algorithmus, in dem der Vorschlag pdf$q()$ist nicht symmetrisch liegt an Hasting (1970) und aus diesem Grund wird der Algorithmus oft auch Metropolis-Hasting genannt. Darüber hinaus wurde hier der globale Metropolis-Algorithmus im Gegensatz zum lokalen Algorithmus beschrieben, bei dem ein Zyklus nur eine Komponente von beeinflusst${\bf x}$.

Der Wikipedia-Artikel ist eine gute Ergänzung zu lesen. Wie Sie sehen können, hat die Metropolis auch ein "Korrekturverhältnis", aber, wie oben erwähnt, hat Hastings eine Modifikation eingeführt, die nicht symmetrische Angebotsverteilungen ermöglicht.

Der Metropolis-Algorithmus ist im R-Paket mcmcunter dem Befehl implementiert metrop().

Andere Codebeispiele:

http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/metrop/

http://pcl.missouri.edu/jeff/node/322

http://darrenjw.wordpress.com/2010/08/15/metropolis-hastings-mcmc-algorithms/