Comprendre Metropolis-Hastings avec une distribution de proposition asymétrique

J'ai essayé de comprendre l'algorithme de Metropolis-Hastings afin d'écrire un code pour estimer les paramètres d'un modèle (ie $f(x)=a*x$ ). Selon la bibliographie, l'algorithme de Metropolis-Hastings comprend les étapes suivantes:

• Générer $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}$

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

Comment je voudrais poser quelques questions:

• La bibliographie indique que si est une distribution symétrique, le rapport devient 1 et l'algorithme est appelé Metropolis. Est-ce exact? La seule différence entre Metropolis et Metropolis-Hastings est que la première utilise une distribution symétrique? Qu'en est-il de la "Random Walk" Metropolis (-Hastings)? En quoi diffère-t-il des deux autres?$q$$q(x|y)/q(y|x)$
• La plupart des exemples de code que je trouve en ligne utilisent une distribution de proposition gaussienne et donc où est le rapport de vraisemblance. Et si la distribution de la proposition est une distribution de Poisson? Je pense que je comprends rationnellement pourquoi ce rapport ne devient pas 1 lorsque j'utilise une distribution asymétrique, mais je ne suis pas sûr de bien le comprendre mathématiquement ou comment l'implémenter avec du code. Quelqu'un pourrait-il me donner un exemple de code simple (C, python, R, pseudo-code ou tout ce que vous préférez) de l'algorithme Metropolis-Hastings utilisant une distribution de proposition non symétrique?$q$$\rho(x,y)=\min( f(y)/f(x),1)$$f(y)/f(x)$

La bibliographie indique que si q est une distribution symétrique, le rapport q (x | y) / q (y | x) devient 1 et l'algorithme est appelé Metropolis. Est-ce exact?

Oui, c'est correct. L'algorithme Metropolis est un cas particulier de l'algorithme MH.

Qu'en est-il de la "Random Walk" Metropolis (-Hastings)? En quoi diffère-t-il des deux autres?

Dans une marche aléatoire, la distribution de la proposition est recentrée après chaque étape à la dernière valeur générée par la chaîne. Généralement, dans une marche aléatoire, la distribution de la proposition est gaussienne, auquel cas cette marche aléatoire satisfait l'exigence de symétrie et l'algorithme est une métropole. Je suppose que vous pourriez effectuer une marche "pseudo" aléatoire avec une distribution asymétrique qui ferait dériver trop les propositions dans la direction opposée de l'inclinaison (une distribution asymétrique à gauche favoriserait les propositions vers la droite). Je ne sais pas pourquoi vous le feriez, mais vous pourriez et ce serait un algorithme de hastings de métropole (c'est-à-dire nécessiter le terme de ratio supplémentaire).

En quoi diffère-t-il des deux autres?

Dans un algorithme de marche non aléatoire, les distributions de proposition sont fixes. Dans la variante de marche aléatoire, le centre de la distribution de proposition change à chaque itération.

Et si la distribution de la proposition est une distribution de Poisson?

Ensuite, vous devez utiliser MH au lieu de simplement métropole. Vraisemblablement, cela consisterait à échantillonner une distribution discrète, sinon vous ne voudriez pas utiliser une fonction discrète pour générer vos propositions.

Dans tous les cas, si la distribution d'échantillonnage est tronquée ou si vous avez une connaissance préalable de son asymétrie, vous souhaiterez probablement utiliser une distribution d'échantillonnage asymétrique et devrez donc utiliser des métropoles-hastings.

Quelqu'un pourrait-il me donner un exemple de code simple (C, python, R, pseudo-code ou tout ce que vous préférez)?

Voici la métropole:

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

Essayons de l'utiliser pour échantillonner une distribution bimodale. Voyons d'abord ce qui se passe si nous utilisons une marche aléatoire pour notre proposition:

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

Essayons maintenant d'échantillonner en utilisant une distribution de proposition fixe et voyons ce qui se passe:

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

Cela semble correct au premier abord, mais si nous examinons la densité postérieure ...

plot(density(y_trace2))

nous verrons qu'il est complètement bloqué au maximum local. Ce n'est pas tout à fait surprenant puisque nous avons en fait centré notre distribution de propositions là-bas. La même chose se produit si nous centrons cela sur l'autre mode:

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

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

Nous pouvons essayer de laisser tomber notre proposition entre les deux modes, mais nous devrons régler la variance très haut pour avoir une chance d'explorer l'un d'eux

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! 

Remarquez comment le choix du centre de distribution de notre proposition a un impact significatif sur le taux d'acceptation de notre échantillonneur.

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

Nous encore coincés dans le plus proche des deux modes. Essayons de le déposer directement entre les deux modes.

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

Enfin, nous nous rapprochons de ce que nous recherchions. Théoriquement, si nous laissons l'échantillonneur fonctionner suffisamment longtemps, nous pouvons obtenir un échantillon représentatif de l'une de ces distributions de proposition, mais la marche aléatoire a produit un échantillon utilisable très rapidement, et nous avons dû tirer parti de notre connaissance de la façon dont le postérieur était supposé chercher à régler les distributions d'échantillonnage fixes pour produire un résultat utilisable (ce qui, à vrai dire, nous ne l'avons pas encore tout à fait y_trace4).

Je vais essayer de mettre à jour avec un exemple de hastings métropole plus tard. Vous devriez être en mesure de voir assez facilement comment modifier le code ci-dessus pour produire un algorithme de hastings de métropole (indice: il vous suffit d'ajouter le rapport supplémentaire dans le logRcalcul).