Путать с вариациями MCMC Метрополис-Гастингс: Случайное прохождение, Неслучайное прохождение, Независимое, Метрополис

15

За последние несколько недель я пытался понять MCMC и алгоритм (ы) Метрополис-Гастингс. Каждый раз, когда я думаю, что понимаю это, я понимаю, что я неправ. Большинство примеров кода, которые я нахожу онлайн, реализуют то, что не согласуется с описанием. то есть: они говорят, что они реализуют Метрополис-Гастингс, но на самом деле они реализуют мегаполис с произвольной ходьбой. Другие (почти всегда) молча пропускают реализацию поправочного коэффициента Гастингса, потому что они используют симметричное распределение предложений. На самом деле, я не нашел ни одного простого примера, который бы рассчитывал соотношение. Это делает меня еще более запутанным. Может ли кто-нибудь дать мне примеры кода (на любом языке) следующего:

  • Алгоритм Метрополиса-Гастингса с не случайным блужданием по Ваниле с вычислением поправочного коэффициента Гастингса (даже если при использовании симметричного распределения предложений это будет 1).
  • Ванильный алгоритм случайного блуждания Метрополиса-Гастингса.
  • Ванильный алгоритм Метрополиса-Гастингса.

Не нужно предоставлять алгоритмы Метрополиса, потому что, если я не ошибаюсь, единственное различие между Метрополисом и Метрополисом-Гастингсом состоит в том, что первые всегда выбирают из симметричного распределения и, следовательно, у них нет поправочного коэффициента Гастингса. Не нужно давать подробное объяснение алгоритмов. Я понимаю основы, но меня смущают разные названия для разных вариаций алгоритма Метрополиса-Гастингса, а также то, как вы практически реализуете поправочный коэффициент Гастингса на MH не случайного блуждания Vanilla. Пожалуйста, не копируйте вставки ссылок, которые частично отвечают на мои вопросы, потому что, скорее всего, я их уже видел. Эти ссылки привели меня в замешательство. Спасибо.

эстрона
источник

Ответы:

10

Вот, пожалуйста, три примера. Я сделал код гораздо менее эффективным, чем в реальном приложении, чтобы сделать логику более понятной (надеюсь).

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

Для ванильного пробоотборника мы получаем:

> 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 

что является малой вероятностью принятия, но все же ... настройка предложения поможет здесь, или принятие другого. Вот результаты предложения случайной прогулки:

> 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 

Подобные результаты, как можно было бы надеяться, и лучшая вероятность принятия (нацеливание на ~ 50% с одним параметром.)

И, для полноты, сэмплер независимости:

> 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 

Поскольку он не «приспосабливается» к форме задней части тела, он, как правило, имеет наименьшую вероятность принятия и его сложнее всего настроить для этой проблемы.

Обратите внимание, что, вообще говоря, мы бы предпочли предложения с более толстыми хвостами, но это совсем другая тема.

jbowman
источник
Q
1
@floyd - это полезно в ряде ситуаций, например, если у вас есть приличное представление о местонахождении центра распределения (например, потому что вы рассчитали оценки MLE или MOM) и можете выбрать предложение с полным хвостом распределение, или если время вычислений на итерацию очень мало, и в этом случае вы можете запустить очень длинную цепочку (что компенсирует низкие показатели приемлемости) - таким образом, вы сэкономите время на анализ и программирование, которое может быть намного больше, чем даже неэффективное время выполнения. Это не было бы типичным предложением с первой попытки, хотя, скорее всего, это была бы случайная прогулка.
jbowman
Qп(ИксT+1|ИксT)
1
п(ИксT+1|ИксT)знак равноп(ИксT+1)
1

Видеть:

Q()Икс

Статья Википедии является хорошим комплементарным чтением. Как вы можете видеть, у Метрополиса также есть «коэффициент коррекции», но, как упоминалось выше, Гастингс представил модификацию, которая допускает несимметричное распределение предложений.

Алгоритм Метрополис реализован в пакете R mcmcпод командой metrop().

Другие примеры кода:

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/

Фриц Ланг
источник
Спасибо за ваш ответ. К сожалению, он не отвечает ни на один из моих вопросов. Я вижу только мегаполис с произвольной ходьбой, мегаполис с неслучайной прогулкой и независимый МЗ Коэффициент коррекции Гастингса dnorm(can,mu,sig)/dnorm(x,mu,sig)в сэмплере независимости первой ссылки не равен 1. Я думал, что он должен был быть равен 1 при использовании симметричного распределения предложений. Это потому, что это независимый пробоотборник, а не обычный MH не-случайного блуждания? Если да, то каков коэффициент Гастингса для простого MH без случайного блуждания?
AstrOne
п(текущий|предложение)знак равноп(предложение|текущий) , т. Е. Если отношение не равно 1, вы не можете их игнорировать.
Jbowman