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