Mylić z odmianami MCMC Metropolis-Hastings: Random-Walk, Non-Random-Walk, Independent, Metropolis

15

W ciągu ostatnich kilku tygodni starałem się zrozumieć MCMC i algorytm (y) Metropolis-Hastings. Za każdym razem, gdy myślę, że to rozumiem, zdaję sobie sprawę, że się mylę. Większość przykładów kodu, które znajduję w Internecie, implementuje coś, co nie jest zgodne z opisem. tj.: Mówią, że wdrażają Metropolis-Hastings, ale w rzeczywistości wdrażają metropolię o swobodnym marszu. Inni (prawie zawsze) po cichu pomijają implementację współczynnika korekcji Hastingsa, ponieważ używają symetrycznego rozkładu propozycji. Właściwie nie znalazłem ani jednego prostego przykładu, który oblicza stosunek do tej pory. To sprawia, że ​​jestem jeszcze bardziej zdezorientowany. Czy ktoś może mi podać przykłady kodu (w dowolnym języku):

  • Waniliowy algorytm obliczania współczynnika korekcji Hastingsa metodą losowego spaceru metropolii-Hastingsa (nawet jeśli przy zastosowaniu symetrycznego rozkładu propozycji będzie to 1)
  • Algorytm Vanilla Random Walk Metropolis-Hastings.
  • Algorytm Vanilla Independent Metropolis-Hastings.

Nie trzeba podawać algorytmów Metropolis, ponieważ jeśli się nie mylę, jedyną różnicą między Metropolis i Metropolis-Hastings jest to, że pierwsze z nich zawsze próbkują z rozkładu symetrycznego, a zatem nie mają współczynnika korekcji Hastingsa. Nie trzeba szczegółowo wyjaśniać algorytmów. Rozumiem podstawy, ale jestem trochę mylony ze wszystkimi różnymi nazwami dla różnych odmian algorytmu Metropolis-Hastings, ale także z tym, jak praktycznie implementujesz współczynnik korekcji Hastingsa w waniliowym MH nieprzypadkowym. Nie kopiuj wklejonych linków, które częściowo odpowiadają na moje pytania, ponieważ najprawdopodobniej już je widziałem. Te linki doprowadziły mnie do tego zamieszania. Dziękuję Ci.

AstrOne
źródło

Odpowiedzi:

10

Proszę bardzo - trzy przykłady. Uczyniłem kod o wiele mniej wydajnym niż w prawdziwej aplikacji, aby logika była bardziej przejrzysta (mam nadzieję).

# We'll assume estimation of a Poisson mean as a function of x
x <- runif(100)
y <- rpois(100,5*x)  # beta = 5 where mean(y[i]) = beta*x[i]

# Prior distribution on log(beta): t(5) with mean 2 
# (Very spread out on original scale; median = 7.4, roughly)
log_prior <- function(log_beta) dt(log_beta-2, 5, log=TRUE)

# Log likelihood
log_lik <- function(log_beta, y, x) sum(dpois(y, exp(log_beta)*x, log=TRUE))

# Random Walk Metropolis-Hastings 
# Proposal is centered at the current value of the parameter

rw_proposal <- function(current) rnorm(1, current, 0.25)
rw_p_proposal_given_current <- function(proposal, current) dnorm(proposal, current, 0.25, log=TRUE)
rw_p_current_given_proposal <- function(current, proposal) dnorm(current, proposal, 0.25, log=TRUE)

rw_alpha <- function(proposal, current) {
   # Due to the structure of the rw proposal distribution, the rw_p_proposal_given_current and
   # rw_p_current_given_proposal terms cancel out, so we don't need to include them - although
   # logically they are still there:  p(prop|curr) = p(curr|prop) for all curr, prop
   exp(log_lik(proposal, y, x) + log_prior(proposal) - log_lik(current, y, x) - log_prior(current))
}

# Independent Metropolis-Hastings
# Note: the proposal is independent of the current value (hence the name), but I maintain the
# parameterization of the functions anyway.  The proposal is not ignorable any more
# when calculation the acceptance probability, as p(curr|prop) != p(prop|curr) in general.

ind_proposal <- function(current) rnorm(1, 2, 1) 
ind_p_proposal_given_current <- function(proposal, current) dnorm(proposal, 2, 1, log=TRUE)
ind_p_current_given_proposal <- function(current, proposal) dnorm(current, 2, 1, log=TRUE)

ind_alpha <- function(proposal, current) {
   exp(log_lik(proposal, y, x)  + log_prior(proposal) + ind_p_current_given_proposal(current, proposal) 
       - log_lik(current, y, x) - log_prior(current) - ind_p_proposal_given_current(proposal, current))
}

# Vanilla Metropolis-Hastings - the independence sampler would do here, but I'll add something
# else for the proposal distribution; a Normal(current, 0.1+abs(current)/5) - symmetric but with a different
# scale depending upon location, so can't ignore the proposal distribution when calculating alpha as
# p(prop|curr) != p(curr|prop) in general

van_proposal <- function(current) rnorm(1, current, 0.1+abs(current)/5)
van_p_proposal_given_current <- function(proposal, current) dnorm(proposal, current, 0.1+abs(current)/5, log=TRUE)
van_p_current_given_proposal <- function(current, proposal) dnorm(current, proposal, 0.1+abs(proposal)/5, log=TRUE)

van_alpha <- function(proposal, current) {
   exp(log_lik(proposal, y, x)  + log_prior(proposal) + ind_p_current_given_proposal(current, proposal) 
       - log_lik(current, y, x) - log_prior(current) - ind_p_proposal_given_current(proposal, current))
}


# Generate the chain
values <- rep(0, 10000) 
u <- runif(length(values))
naccept <- 0
current <- 1  # Initial value
propfunc <- van_proposal  # Substitute ind_proposal or rw_proposal here
alphafunc <- van_alpha    # Substitute ind_alpha or rw_alpha here
for (i in 1:length(values)) {
   proposal <- propfunc(current)
   alpha <- alphafunc(proposal, current)
   if (u[i] < alpha) {
      values[i] <- exp(proposal)
      current <- proposal
      naccept <- naccept + 1
   } else {
      values[i] <- exp(current)
   }
}
naccept / length(values)
summary(values)

W przypadku waniliowego samplera otrzymujemy:

> naccept / length(values)
[1] 0.1737
> summary(values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.843   5.153   5.388   5.378   5.594   6.628 

co jest niskim prawdopodobieństwem akceptacji, ale nadal ... dostrojenie propozycji pomogłoby tutaj lub przyjęcie innej. Oto wyniki propozycji losowego spaceru:

> naccept / length(values)
[1] 0.2902
> summary(values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.718   5.147   5.369   5.370   5.584   6.781 

Podobne wyniki, jak można się spodziewać, i większe prawdopodobieństwo akceptacji (dążenie do ~ 50% z jednym parametrem).

I, dla kompletności, próbnik niezależności:

> naccept / length(values)
[1] 0.0684
> summary(values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  3.990   5.162   5.391   5.380   5.577   8.802 

Ponieważ nie „dostosowuje się” do kształtu tylnej części ciała, ma najgorsze prawdopodobieństwo akceptacji i najtrudniej jest dostroić się do tego problemu.

Zauważ, że ogólnie rzecz biorąc wolelibyśmy propozycje z grubszymi ogonami, ale to zupełnie inny temat.

łucznik
źródło
Q
1
@floyd - jest przydatny w wielu sytuacjach, na przykład, jeśli masz dobre wyobrażenie o lokalizacji centrum dystrybucji (np. ponieważ obliczyłeś szacunki MLE lub MOM) i możesz wybrać grubą propozycję dystrybucja lub jeśli czas obliczeń na iterację jest bardzo krótki, w takim przypadku możesz uruchomić bardzo długi łańcuch (który rekompensuje niskie wskaźniki akceptacji) - oszczędzając w ten sposób czas analizy i programowania, który może być znacznie dłuższy niż nawet nieefektywne środowisko uruchomieniowe. Nie byłaby to jednak typowa propozycja pierwszej próby, prawdopodobnie byłby to przypadkowy spacer.
jbowman
Qp(xt+1|xt)
1
p(xt+1|xt)=p(xt+1)
1

Widzieć:

q()x

Artykuł w Wikipedii jest dobrym uzupełnieniem. Jak widać, Metropolis ma również „współczynnik korekcji”, ale, jak wspomniano powyżej, Hastings wprowadził modyfikację, która pozwala na niesymetryczne rozkłady propozycji.

Algorytm Metropolis zaimplementowana w pakiecie R mcmcna podstawie polecenia metrop().

Inne przykłady kodu:

http://www.mas.ncl.ac.uk/~ndjw1/teaching/sim/metrop/

http://pcl.missouri.edu/jeff/node/322

http://darrenjw.wordpress.com/2010/08/15/metropolis-hastings-mcmc-algorithms/

Fritz Lang
źródło
Dziękuję za odpowiedź. Niestety nie odpowiada na żadne z moich pytań. Widzę tylko metropolię o losowym marszu, metropolię bez losowego chodu i niezależne MH. Współczynnik korekcji Hastingsa dnorm(can,mu,sig)/dnorm(x,mu,sig)w próbniku niezależności pierwszego łącza nie jest równy 1. Pomyślałem, że powinien on być równy 1, gdy używa się symetrycznego rozkładu propozycji. Czy to dlatego, że jest to niezależny sampler, a nie zwykły MH nieprzypadkowy? Jeśli tak, jaki jest stosunek Hastingsa dla zwykłego MH bez marszu?
AstrOne,
p(obecny|wniosek)=p(wniosek|obecny)