Понимание Metropolis-Hastings с асимметричным распределением предложений

Я пытался понять алгоритм Метрополиса-Гастингса, чтобы написать код для оценки параметров модели (т.е. $f(x)=a*x$ ). Согласно библиографии алгоритм Метрополиса-Гастингса имеет следующие этапы:

• Сгенерировать $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)$

Как я хотел бы задать несколько вопросов:

• Библиография утверждает, что если - симметричное распределение, отношение q ( x | y ) / q ( y | x ) становится равным 1, и алгоритм называется Metropolis. Это верно? Единственная разница между Метрополисом и Метрополисом-Гастингсом заключается в том, что в первом используется симметричное распределение? А как насчет "Случайной прогулки" Метрополис (-Хастингс)? Чем он отличается от двух других?$q$$q(x|y)/q(y|x)$
• Большая часть примера кода, который я нахожу в Интернете, использует гауссово распределение предложений и, таким образом, ρ ( x , y ) = min ( f ( y ) / f ( x ) , 1 ) где f ( y ) / f ( x $q$$\rho(x,y)=\min( f(y)/f(x),1)$$f(y)/f(x)$это отношение правдоподобия. Что если распределение предложений является распределением Пуассона? Я думаю, что понимаю рационально, почему это соотношение не становится равным 1 при использовании асимметричного распределения, но я не совсем уверен, понимаем ли это математически или как реализовать его с помощью кода. Может ли кто-нибудь дать мне простой пример кода (C, python, R, псевдокод или что-то другое) алгоритма Метрополиса-Гастингса, использующего несимметричное распределение предложений?

Библиография утверждает, что если q - симметричное распределение, отношение q (x | y) / q (y | x) становится равным 1, и алгоритм называется Metropolis. Это верно?

Да, это правильно. Алгоритм Метрополиса является частным случаем алгоритма MH.

А как насчет "Случайной прогулки" Метрополис (-Хастингс)? Чем он отличается от двух других?

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

Чем он отличается от двух других?

В алгоритме неслучайного блуждания распределения предложений фиксированы. В варианте случайного обхода центр распределения предложения меняется на каждой итерации.

Что если распределение предложений является распределением Пуассона?

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

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

Может ли кто-нибудь дать мне простой пример кода (C, Python, R, псевдокод или как вы предпочитаете)?

Вот мегаполис:

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

Давайте попробуем использовать это для демонстрации бимодального распределения. Во-первых, давайте посмотрим, что произойдет, если мы используем случайную прогулку для нашего propsal:

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

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

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

Сначала это выглядит нормально, но если мы посмотрим на заднюю плотность ...

plot(density(y_trace2))

мы увидим, что он полностью застрял на локальном максимуме. Это не совсем удивительно, так как мы фактически сосредоточили свое распространение там. То же самое происходит, если мы центрируем это в другом режиме:

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

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

Мы можем попытаться отбросить наше предложение между двумя режимами, но нам нужно будет установить очень высокую дисперсию, чтобы иметь возможность исследовать любой из них.

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! 

Обратите внимание, как выбор центра распространения нашего предложения оказывает существенное влияние на скорость принятия нашего пробоотборника.

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

Мы все еще застряли в ближе двух режимов. Давайте попробуем сбросить это прямо между двумя режимами.

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

Наконец, мы приближаемся к тому, что искали. Теоретически, если мы позволим сэмплеру работать достаточно долго, мы сможем получить репрезентативную выборку из любого из этих распределений предложений, но случайное блуждание произвело пригодную для использования выборку очень быстро, и нам пришлось воспользоваться нашими знаниями о том, как предполагался апостериорный попытаться настроить распределения фиксированной выборки, чтобы получить пригодный для использования результат (который, честно говоря, у нас пока нет y_trace4).

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