Memahami Metropolis-Hastings dengan distribusi proposal asimetris

Saya telah mencoba memahami algoritma Metropolis-Hastings untuk menulis kode untuk memperkirakan parameter model (yaitu ). Menurut bibliografi, algoritma Metropolis-Hastings memiliki langkah-langkah berikut:$f(x)=a*x$

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

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

Bagaimana saya ingin mengajukan beberapa pertanyaan:

• Daftar pustaka menyatakan bahwa jika $q$ adalah distribusi simetris rasio $q(x|y)/q(y|x)$ menjadi 1 dan algoritma ini disebut Metropolis. Apakah itu benar? Satu-satunya perbedaan antara Metropolis dan Metropolis-Hastings adalah yang pertama menggunakan distribusi simetris? Bagaimana dengan Metropolis "Random Walk" (-Hastings)? Apa bedanya dengan dua lainnya?
• Sebagian besar kode contoh yang saya temukan on-line menggunakan distribusi proposal Gaussian dan karenanya mana adalah rasio kemungkinan. Bagaimana jika distribusi proposal adalah distribusi Poisson? Saya pikir saya mengerti secara rasional mengapa rasio itu tidak menjadi 1 ketika menggunakan distribusi asimetris tetapi saya tidak yakin apakah memahaminya secara matematis atau bagaimana menerapkannya dengan kode. Bisakah seseorang memberi saya kode sederhana (C, python, R, pseudo-code atau apa pun yang Anda suka) contoh dari algoritma Metropolis-Hastings menggunakan distribusi proposal yang tidak simetris?$q$$\rho(x,y)=\min( f(y)/f(x),1)$$f(y)/f(x)$

Daftar pustaka menyatakan bahwa jika q adalah distribusi simetris rasio q (x | y) / q (y | x) menjadi 1 dan algoritma ini disebut Metropolis. Apakah itu benar?

Ya ini benar. Algoritma Metropolis adalah kasus khusus dari algoritma MH.

Bagaimana dengan Metropolis "Random Walk" (-Hastings)? Apa bedanya dengan dua lainnya?

Secara acak, distribusi proposal dipusatkan kembali setelah setiap langkah pada nilai terakhir yang dihasilkan oleh rantai. Secara umum, dalam jalan acak distribusi proposal adalah gaussian, dalam hal ini jalan acak memenuhi persyaratan simetri dan algoritmanya adalah metropolis. Saya kira Anda dapat melakukan jalan acak "semu" dengan distribusi asimetris yang akan menyebabkan proposal terlalu melayang ke arah yang berlawanan dari kemiringan (distribusi condong ke kiri akan mendukung proposal ke kanan). Saya tidak yakin mengapa Anda melakukan ini, tetapi Anda bisa dan itu akan menjadi algoritma metropolis hastings (yaitu memerlukan istilah rasio tambahan).

Apa bedanya dengan dua lainnya?

Dalam algoritma berjalan non-acak, distribusi proposal diperbaiki. Dalam varian random walk, pusat distribusi proposal berubah pada setiap iterasi.

Bagaimana jika distribusi proposal adalah distribusi Poisson?

Maka Anda perlu menggunakan MH bukan hanya metropolis. Agaknya ini akan menjadi sampel distribusi diskrit, jika tidak, Anda tidak akan ingin menggunakan fungsi diskrit untuk menghasilkan proposal Anda.

Dalam hal apa pun, jika distribusi sampel dipotong atau Anda memiliki pengetahuan sebelumnya tentang kemiringannya, Anda mungkin ingin menggunakan distribusi sampling asimetris dan karenanya perlu menggunakan metropolis-hastings.

Bisakah seseorang memberi saya kode sederhana (C, python, R, pseudo-code atau apa pun yang Anda suka) misalnya?

Inilah kota metropolis:

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

Mari kita coba gunakan ini untuk mencicipi distribusi bimodal. Pertama, mari kita lihat apa yang terjadi jika kita menggunakan jalan acak untuk propsal kita:

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

Sekarang mari kita coba pengambilan sampel menggunakan distribusi proposal tetap dan lihat apa yang terjadi:

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

Ini terlihat ok pada awalnya, tetapi jika kita melihat kepadatan posterior ...

plot(density(y_trace2))

kita akan melihat bahwa itu benar-benar macet maksimum lokal. Ini tidak sepenuhnya mengejutkan karena kami sebenarnya memusatkan distribusi proposal kami di sana. Hal yang sama terjadi jika kita memusatkan ini pada mode lain:

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

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

Kita dapat mencoba membuang proposal kita di antara dua mode, tetapi kita perlu mengatur varians sangat tinggi untuk memiliki kesempatan menjelajahi salah satu dari mereka.

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! 

Perhatikan bagaimana pilihan pusat distribusi proposal kami memiliki dampak signifikan pada tingkat penerimaan sampler kami.

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

Kami masih terjebak di dekat dua mode. Mari kita coba menjatuhkan ini langsung antara dua mode.

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

Akhirnya, kami semakin dekat dengan apa yang kami cari. Secara teoritis, jika kita membiarkan sampler bekerja cukup lama kita bisa mendapatkan sampel yang representatif dari salah satu dari distribusi proposal ini, tetapi jalan acak menghasilkan sampel yang dapat digunakan dengan sangat cepat, dan kita harus mengambil keuntungan dari pengetahuan kita tentang bagaimana posterior itu seharusnya untuk melihat untuk menyetel distribusi sampling tetap untuk menghasilkan hasil yang dapat digunakan (yang, sejujurnya, kami belum cukup dalam y_trace4).

Saya akan mencoba memperbarui dengan contoh metropolis hastings nanti. Anda harus dapat melihat dengan cukup mudah bagaimana memodifikasi kode di atas untuk menghasilkan algoritma metropolis hastings (petunjuk: Anda hanya perlu menambahkan rasio tambahan ke dalam logRperhitungan).