Compreendendo Metropolis-Hastings com distribuição de proposta assimétrica

Eu tenho tentado entender o algoritmo Metropolis-Hastings para escrever um código para estimar os parâmetros de um modelo (ou seja, ). Segundo a bibliografia, o algoritmo Metropolis-Hastings possui as seguintes etapas:$f(x)=a*x$

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

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

Como eu gostaria de fazer algumas perguntas:

• A bibliografia afirma que, se é uma distribuição simétrica, a razão q ( x | y ) / q ( y | x ) se torna 1 e o algoritmo é chamado Metropolis. Isso está correto? A única diferença entre Metropolis e Metropolis-Hastings é que o primeiro usa uma distribuição simétrica? E quanto à metrópole "Random Walk" (-Hastings)? Como isso difere dos outros dois?$q$$q(x|y)/q(y|x)$
• A maioria do código de exemplo que eu acho on-line usa uma distribuição de proposta gaussiana e, portanto, ρ ( x , y ) = min ( f ( y ) / f ( x ) , 1 ) onde f ( y ) / f ( x $q$$\rho(x,y)=\min( f(y)/f(x),1)$$f(y)/f(x)$é a razão de verossimilhança. E se a distribuição da proposta for uma distribuição de Poisson? Acho que entendo racionalmente por que essa proporção não se torna 1 ao usar uma distribuição assimétrica, mas não tenho certeza se a entendo matematicamente ou como implementá-la com código. Alguém poderia me dar um exemplo simples de código (C, python, R, pseudo-código ou o que você preferir) do algoritmo Metropolis-Hastings usando uma distribuição de proposta não simétrica?

A bibliografia afirma que, se q é uma distribuição simétrica, a razão q (x | y) / q (y | x) se torna 1 e o algoritmo é chamado Metropolis. Isso está correto?

Sim isto está correcto. O algoritmo Metropolis é um caso especial do algoritmo MH.

E quanto à metrópole "Random Walk" (-Hastings)? Como isso difere dos outros dois?

Em uma caminhada aleatória, a distribuição da proposta é centrada novamente após cada etapa no último valor gerado pela cadeia. Geralmente, em uma caminhada aleatória, a distribuição da proposta é gaussiana; nesse caso, essa caminhada aleatória atende ao requisito de simetria e o algoritmo é metrópole. Suponho que você possa executar um passeio "pseudo" aleatório com uma distribuição assimétrica, o que faria com que as propostas se desviassem demais na direção oposta à inclinação (uma distribuição inclinada à esquerda favoreceria propostas à direita). Não sei por que você faria isso, mas você poderia e seria um algoritmo de detenções de metrópoles (por exemplo, exigir o termo de proporção adicional).

Como isso difere dos outros dois?

Em um algoritmo de caminhada não aleatória, as distribuições da proposta são fixas. Na variante de passeio aleatório, o centro da distribuição da proposta é alterado a cada iteração.

E se a distribuição da proposta for uma distribuição de Poisson?

Então você precisa usar MH em vez de apenas metrópole. Presumivelmente, isso seria uma amostra de uma distribuição discreta; caso contrário, você não desejaria usar uma função discreta para gerar suas propostas.

De qualquer forma, se a distribuição de amostragem estiver truncada ou você tiver conhecimento prévio de sua inclinação, provavelmente desejará usar uma distribuição de amostragem assimétrica e, portanto, precisará usar hastings de metrópoles.

Alguém poderia me dar um exemplo simples de código (C, python, R, pseudo-código ou o que você preferir)?

Aqui está a metrópole:

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

Vamos tentar usar isso para provar uma distribuição bimodal. Primeiro, vamos ver o que acontece se usarmos uma caminhada aleatória para nossa 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

Agora, vamos tentar fazer uma amostragem usando uma distribuição de proposta fixa e ver o que acontece:

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

Parece bom no começo, mas se dermos uma olhada na densidade posterior ...

plot(density(y_trace2))

veremos que está completamente preso no máximo local. Isso não é totalmente surpreendente, pois na verdade centralizamos nossa distribuição de propostas lá. O mesmo acontece se centralizarmos isso no outro modo:

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

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

Podemos tentar abandonar nossa proposta entre os dois modos, mas precisaremos definir uma variação muito alta para ter a chance de explorar qualquer um deles

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 como a escolha do centro da distribuição da proposta tem um impacto significativo na taxa de aceitação do nosso amostrador.

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

Nós ainda ficar preso na mais perto dos dois modos. Vamos tentar deixar isso diretamente entre os dois 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, estamos nos aproximando do que estávamos procurando. Teoricamente, se deixarmos o amostrador funcionar o tempo suficiente, podemos obter uma amostra representativa de qualquer uma dessas distribuições de propostas, mas a caminhada aleatória produziu uma amostra utilizável muito rapidamente e tivemos que tirar proveito de nosso conhecimento de como o posterior era suposto. procurar ajustar as distribuições de amostragem fixas para produzir um resultado utilizável (que, verdade seja dita, ainda não temos y_trace4).

Vou tentar atualizar com um exemplo de pressa nas metrópoles mais tarde. Você deve ver com bastante facilidade como modificar o código acima para produzir um algoritmo de hastings de metrópoles (dica: você só precisa adicionar a taxa suplementar ao logRcálculo).