Zrozumienie Metropolis-Hastings z asymetrycznym rozkładem propozycji

Próbowałem zrozumieć algorytm Metropolis-Hastings, aby napisać kod do oszacowania parametrów modelu (tj. ). Według bibliografii algorytm Metropolis-Hastings ma następujące kroki:$f(x)=a*x$

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

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

Jak chciałbym zadać kilka pytań:

• Bibliografia stwierdza, że ​​jeśli $q$ jest rozkładem symetrycznym, stosunek $q(x|y)/q(y|x)$ wynosi 1, a algorytm nazywa się Metropolis. Czy to jest poprawne? Jedyną różnicą między Metropolis i Metropolis-Hastings jest to, że w pierwszej zastosowano rozkład symetryczny? A co z Metropolis „Random Walk” (-Hastings)? Czym różni się od pozostałych dwóch?
• Większość przykładowego kodu, który znajduję on-line, wykorzystuje rozkład propozycji Gaussa $q$ a zatem $\rho(x,y)=\min( f(y)/f(x),1)$ gdzie $f(y)/f(x)$ to współczynnik prawdopodobieństwa. Co jeśli rozkład propozycji jest rozkładem Poissona? Wydaje mi się, że rozumiem racjonalnie, dlaczego ten stosunek nie staje się 1, gdy stosuje się rozkład asymetryczny, ale nie jestem pewien, czy rozumiem go matematycznie lub jak zaimplementować go za pomocą kodu. Czy ktoś mógłby mi podać prosty kod (C, python, R, pseudo-kod lub cokolwiek innego) przykład algorytmu Metropolis-Hastings wykorzystujący niesymetryczny rozkład propozycji?

Bibliografia stwierdza, że ​​jeśli q jest rozkładem symetrycznym, stosunek q (x | y) / q (y | x) wynosi 1, a algorytm nazywa się Metropolis. Czy to jest poprawne?

Tak, to jest poprawne. Algorytm Metropolis jest szczególnym przypadkiem algorytmu MH.

A co z Metropolis „Random Walk” (-Hastings)? Czym różni się od pozostałych dwóch?

W przypadkowym spacerze rozkład propozycji jest ponownie wyśrodkowany po każdym kroku na wartości ostatnio wygenerowanej przez łańcuch. Zasadniczo w przypadku chodzenia losowego rozkład propozycji jest gaussowski, w którym to przypadku chodzenie losowe spełnia wymóg symetrii, a algorytmem jest metropolia. Przypuszczam, że mógłbyś wykonać losowy „pseudo” spacer z asymetrycznym rozkładem, który spowodowałby, że propozycje zbyt dryfowałyby w przeciwnym kierunku przekosu (lewy rozkład pochyliłby propozycje w prawo). Nie jestem pewien, dlaczego to zrobiłeś, ale mógłbyś i byłby to algorytm pospieszny metropolii (tj. Wymagałby dodatkowego współczynnika proporcji).

Czym różni się od pozostałych dwóch?

W przypadku niepolotowego algorytmu chodzenia rozkłady propozycji są ustalone. W wariancie chodzenia losowego środek rozkładu propozycji zmienia się przy każdej iteracji.

Co jeśli rozkład propozycji jest rozkładem Poissona?

Następnie musisz użyć MH zamiast zwykłej metropolii. Przypuszczalnie byłoby to próbkowanie dystrybucji dyskretnej, w przeciwnym razie nie chcesz używać funkcji dyskretnej do generowania swoich propozycji.

W każdym razie, jeśli rozkład próbkowania jest obcięty lub masz wcześniejszą wiedzę na temat jego pochylenia, prawdopodobnie chcesz zastosować asymetryczny rozkład próbkowania, a zatem musisz użyć metropolii.

Czy ktoś mógłby mi podać prosty kod (C, python, R, pseudokod lub cokolwiek innego) przykład?

Oto metropolia:

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

Spróbujmy użyć tego do próbkowania rozkładu bimodalnego. Po pierwsze, zobaczmy, co się stanie, jeśli użyjemy przypadkowego spaceru do naszego rekwizytu:

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

Teraz spróbujmy próbkować przy użyciu ustalonego rozkładu propozycji i zobaczmy, co się stanie:

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

Na początku wygląda to dobrze, ale jeśli spojrzymy na gęstość tylnej ...

plot(density(y_trace2))

zobaczymy, że jest całkowicie zablokowane na lokalnym maksimum. Nie jest to całkowicie zaskakujące, ponieważ właściwie skoncentrowaliśmy tam naszą dystrybucję propozycji. To samo stanie się, jeśli skoncentrujemy to na innym trybie:

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

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

Możemy spróbować upuścić naszą propozycję między dwoma trybami, ale musimy ustawić naprawdę dużą wariancję, aby mieć szansę na zbadanie jednego z nich

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! 

Zauważ, że wybór centrum dystrybucji naszego wniosku ma znaczący wpływ na wskaźnik akceptacji naszego próbnika.

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

My nadal utknąć w bliżej z dwóch trybów. Spróbujmy upuścić to bezpośrednio między dwoma trybami.

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

Wreszcie zbliżamy się do tego, czego szukaliśmy. Teoretycznie, jeśli pozwolimy, aby próbnik działał wystarczająco długo, możemy pobrać reprezentatywną próbkę z dowolnego z tych rozkładów propozycji, ale losowy spacer bardzo szybko wytworzył użyteczną próbkę i musieliśmy skorzystać z naszej wiedzy o tym, jak przypuszczano, że tył aby dostroić ustalone rozkłady próbkowania w celu uzyskania użytecznego wyniku (który, prawdę mówiąc, jeszcze nie mamy y_trace4).

Później postaram się zaktualizować o przykład przyspieszenia metropolii. Powinieneś być w stanie dość łatwo zobaczyć, jak zmodyfikować powyższy kod, aby uzyskać algorytm przyspieszenia metropolii (wskazówka: wystarczy dodać dodatkowy współczynnik do logRobliczeń).