Comprensión de Metropolis-Hastings con distribución de propuesta asimétrica

14

He estado tratando de entender el algoritmo de Metropolis-Hastings para escribir un código para estimar los parámetros de un modelo (es decir, F(X)=unX ). Según la bibliografía, el algoritmo Metropolis-Hastings tiene los siguientes pasos:

  • Generar Ytq(yEl |Xt)
  • Xt+1={Yt,con probabilidadρ(Xt,Yt),Xt,con probabilidad1-ρ(Xt,Yt),

donde ρ(X,y)=min(F(y)F(X)q(XEl |y)q(yEl |X),1)

Cómo me gustaría hacer algunas preguntas:

  • La bibliografía establece que si es una distribución simétrica, la relación q ( x | y ) / q ( y | x ) se convierte en 1 y el algoritmo se llama Metrópolis. ¿Es eso correcto? La única diferencia entre Metropolis y Metropolis-Hastings es que el primero usa una distribución simétrica. ¿Qué pasa con Metrópolis "Random Walk" (-Hastings)? ¿Cómo difiere de los otros dos?qq(XEl |y)/ /q(yEl |X)
  • qρ(X,y)=min(F(y)/ /F(X),1)F(y)/ /F(X)es la razón de probabilidad ¿Qué pasa si la distribución de la propuesta es una distribución de Poisson? Creo que entiendo racionalmente por qué esa relación no se convierte en 1 cuando se usa una distribución asimétrica, pero no estoy muy seguro de si la entiendo matemáticamente o cómo implementarla con código. ¿Podría alguien darme un código simple (C, python, R, pseudocódigo o lo que prefiera) ejemplo del algoritmo Metropolis-Hastings utilizando una distribución de propuesta no simétrica?
AstrOne
fuente
1
Acabo de recordar una excelente publicación de blog sobre un tema relacionado, tal vez esto ayude: darrenjw.wordpress.com/2012/06/04/…
joint_p

Respuestas:

20

La bibliografía establece que si q es una distribución simétrica, la relación q (x | y) / q (y | x) se convierte en 1 y el algoritmo se llama Metrópolis. ¿Es eso correcto?

Si eso es correcto. El algoritmo Metropolis es un caso especial del algoritmo MH.

¿Qué pasa con Metrópolis "Random Walk" (-Hastings)? ¿Cómo difiere de los otros dos?

En una caminata aleatoria, la distribución de la propuesta se vuelve a centrar después de cada paso en el último valor generado por la cadena. Generalmente, en una caminata aleatoria, la distribución de la propuesta es gaussiana, en cuyo caso esta caminata aleatoria satisface el requisito de simetría y el algoritmo es metrópoli. Supongo que podría realizar una caminata "pseudo" aleatoria con una distribución asimétrica que causaría que las propuestas se desvíen demasiado en la dirección opuesta al sesgo (una distribución sesgada a la izquierda favorecería las propuestas hacia la derecha). No estoy seguro de por qué haría esto, pero podría y sería un algoritmo de aceleración de metrópolis (es decir, requiere el término de relación adicional).

¿Cómo difiere de los otros dos?

En un algoritmo de paseo no aleatorio, las distribuciones de la propuesta son fijas. En la variante de caminata aleatoria, el centro de la distribución de la propuesta cambia en cada iteración.

¿Qué pasa si la distribución de la propuesta es una distribución de Poisson?

Entonces necesitas usar MH en lugar de solo metrópolis. Presumiblemente, esto sería probar una distribución discreta, de lo contrario no querría usar una función discreta para generar sus propuestas.

En cualquier caso, si la distribución de muestreo está truncada o si tiene conocimiento previo de su sesgo, es probable que desee utilizar una distribución de muestreo asimétrica y, por lo tanto, debe utilizar hastings de metrópolis.

¿Podría alguien darme un ejemplo de código simple (C, python, R, pseudocódigo o lo que prefiera)?

Aquí está la metrópoli:

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

Intentemos usar esto para probar una distribución bimodal. Primero, veamos qué sucede si usamos una caminata aleatoria para nuestra propsal:

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

ingrese la descripción de la imagen aquí

Ahora intentemos el muestreo usando una distribución de propuesta fija y veamos qué sucede:

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

Esto se ve bien al principio, pero si echamos un vistazo a la densidad posterior ...

plot(density(y_trace2))

ingrese la descripción de la imagen aquí

veremos que está completamente atascado en un máximo local. Esto no es del todo sorprendente, ya que en realidad centramos nuestra distribución de propuestas allí. Lo mismo sucede si centramos esto en el otro modo:

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

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

Podemos intentar dejar nuestra propuesta entre los dos modos, pero necesitaremos establecer la varianza realmente alta para tener la oportunidad de explorar cualquiera de ellos.

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! 

Observe cómo la elección del centro de distribución de nuestra propuesta tiene un impacto significativo en la tasa de aceptación de nuestra muestra.

plot(density(y_trace3))

ingrese la descripción de la imagen aquí

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

Todavía nos quedamos atrapados en el más cercano de los dos modos. Intentemos soltar esto directamente entre los dos modos.

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

ingrese la descripción de la imagen aquí

Finalmente, nos estamos acercando a lo que estábamos buscando. Teóricamente, si dejamos que la muestra funcione lo suficiente, podemos obtener una muestra representativa de cualquiera de estas distribuciones de propuestas, pero la caminata aleatoria produjo una muestra utilizable muy rápidamente, y tuvimos que aprovechar nuestro conocimiento de cómo se suponía que la posterior buscar ajustar las distribuciones de muestreo fijas para producir un resultado utilizable (que, a decir verdad, todavía no tenemos y_trace4).

Intentaré actualizar con un ejemplo de hostigamiento de metrópolis más adelante. Debería poder ver con bastante facilidad cómo modificar el código anterior para producir un algoritmo de aceleración de metrópolis (sugerencia: solo necesita agregar la relación suplementaria en el logRcálculo).

David Marx
fuente
Respuesta impresionante! Muchas gracias! En mi caso, tengo un modelo de parámetros 6-7 y no tengo idea de cómo se vería la distribución posterior (pero podría ser bimodal) porque mis conjuntos de datos a veces son completamente diferentes. Según lo que usted dijo, puedo usar Metropolis (-Hastings) usando una gran variación en la distribución de la propuesta o usar Random Walk Metropolis (-Hastings) con una variación menor en la distribución de la propuesta. Bajo ninguna circunstancia especial, la segunda solución debería converger más rápido a la distribución objetivo. ¿Derecho?
AstrOne
Ahora relacionado con el código Metropolis-Hastings, estaba pensando en reemplazarlo R=exp(logR)con esto: R=exp(logR)*(dnorm(y[t-1],y.prop,my_sigma)/dnorm(y.prop,y[t-1],my_sigma))tanto para la caminata aleatoria como para la no aleatoria MH. ¿Es eso correcto?
AstrOne
1
Básicamente, pero como menciono en el código de la metrópoli: desea realizar sus cálculos en el espacio de registro. Los cálculos de probabilidad tienden a funcionar con valores muy pequeños, por lo que generalmente obtiene resultados mucho mejores agregando logaritmos y exponiendo sus resultados que multiplicando valores brutos juntos.
David Marx
1
yt-1q(ytEl |yt-1)=q(yt)yt-1
1
Usted establece "En un algoritmo de recorrido no aleatorio, las distribuciones de la propuesta son fijas. En la variante de recorrido aleatorio, el centro de la distribución de la propuesta cambia en cada iteración", lo cual no es correcto. Las versiones de MH que no son caminatas aleatorias en su mayoría tienen propuestas que dependen del estado actual de la cadena de Markov, a veces incluso centradas en este estado. Un ejemplo importante es el algoritmo Langevin MCMC. Cuando se corrige la propuesta, este es un algoritmo MH independiente .
Xi'an