Confundido con las variaciones de MCMC Metropolis-Hastings: caminata aleatoria, caminata no aleatoria, independiente, metrópoli

En las últimas semanas, he estado tratando de entender MCMC y los algoritmos de Metropolis-Hastings. Cada vez que creo entenderlo me doy cuenta de que estoy equivocado. La mayoría de los ejemplos de código que encuentro en línea implementan algo que no es consistente con la descripción. es decir: dicen que implementan Metropolis-Hastings pero en realidad implementan metrópolis de paseo aleatorio. Otros (casi siempre) omiten silenciosamente la implementación de la relación de corrección de Hastings porque están utilizando una distribución de propuesta simétrica. En realidad, no he encontrado un solo ejemplo simple que calcule la relación hasta ahora. Eso me confunde aún más. ¿Puede alguien darme ejemplos de código (en cualquier idioma) de lo siguiente:

• Algoritmo de vainilla no aleatorio de Metropolis-Hastings con cálculo de la relación de corrección de Hastings (incluso si esto termina siendo 1 cuando se utiliza una distribución de propuesta simétrica).
• Algoritmo Vanilla Random Walk Metropolis-Hastings.
• Algoritmo Vanilla Independent Metropolis-Hastings.

No es necesario proporcionar los algoritmos de Metropolis porque, si no me equivoco, la única diferencia entre Metropolis y Metropolis-Hastings es que los primeros siempre toman muestras de una distribución simétrica y, por lo tanto, no tienen la relación de corrección de Hastings. No es necesario dar una explicación detallada de los algoritmos. Entiendo los conceptos básicos, pero estoy un poco confundido con todos los diferentes nombres de las diferentes variaciones del algoritmo de Metropolis-Hastings, pero también con la forma en que prácticamente implementa la relación de corrección de Hastings en el MH de caminata no aleatoria de Vanilla. No copie enlaces de pegado que respondan parcialmente a mis preguntas porque lo más probable es que ya los haya visto. Esos enlaces me llevaron a esta confusión. Gracias.

Aquí tienes tres ejemplos. He hecho que el código sea mucho menos eficiente de lo que sería en una aplicación real para aclarar la lógica (espero).

# 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)

Para la muestra de vainilla, obtenemos:

> 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 

lo cual es una baja probabilidad de aceptación, pero aún así ... ajustar la propuesta ayudaría aquí, o adoptar una diferente. Aquí están los resultados de la propuesta de caminata aleatoria:

> 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 

Resultados similares, como uno esperaría, y una mejor probabilidad de aceptación (con el objetivo de ~ 50% con un parámetro).

Y, para completar, la muestra de independencia:

> 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 

Debido a que no se "adapta" a la forma de la parte posterior, tiende a tener la menor probabilidad de aceptación y es más difícil sintonizar bien este problema.

Tenga en cuenta que, en términos generales, preferiríamos propuestas con colas más gruesas, pero ese es un tema completamente diferente.

Ver:

Por construcción, el algoritmo no depende de la constante de normalización, ya que lo que importa es la relación de los pdf. La variación del algoritmo en el que la propuesta pdf$q()$no es simétrico se debe a Hasting (1970) y por esta razón, el algoritmo a menudo también se llama Metropolis-Hasting. Además, lo que se ha descrito aquí es el algoritmo global de Metropolis, en contraste con el local, en el que un ciclo afecta solo a un componente de${\bf x}$.

El artículo de Wikipedia es una buena lectura complementaria. Como puede ver, el Metropolis también tiene una "relación de corrección" pero, como se mencionó anteriormente, Hastings introdujo una modificación que permite distribuciones de propuestas no simétricas.

El algoritmo Metropolis se implementa en el paquete R mcmcbajo el comando metrop().

Otros ejemplos de código:

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/