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

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)=a*x$ ). Según la bibliografía, el algoritmo Metropolis-Hastings tiene los siguientes pasos:

• Generar $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}$

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

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?$q$$q(x|y)/q(y|x)$
• $q$$\rho(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?

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

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

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

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

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