Memahami Metropolis-Hastings dengan distribusi proposal asimetris

14

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)=Sebuahx

  • HasilkanYtq(y|xt)
  • Xt+1={Yt,dengan probabilitasρ(xt,Yt),xt,dengan probabilitas1-ρ(xt,Yt),

di manaρ(x,y)=min(f(y)f(x)q(x|y)q(y|x),1)

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ρ(x,y)=min(f(y)/f(x),1)f(y)/f(x)
AstrOne
sumber
1
Saya baru saja mengingat posting blog yang bagus tentang masalah terkait, mungkin ini membantu: darrenjw.wordpress.com/2012/06/04/…
joint_p

Jawaban:

20

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

masukkan deskripsi gambar di sini

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

masukkan deskripsi gambar di sini

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

masukkan deskripsi gambar di sini

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

masukkan deskripsi gambar di sini

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

David Marx
sumber
Jawaban yang luar biasa! Terima kasih banyak! Dalam kasus saya, saya memiliki model parameter 6-7 dan saya tidak tahu bagaimana distribusi posterior akan terlihat (tetapi mungkin bimodal) karena dataset saya kadang-kadang sangat berbeda. Berdasarkan apa yang Anda katakan, saya bisa menggunakan Metropolis (-Hastings) menggunakan varian yang sangat besar dalam distribusi proposal atau menggunakan Random Walk Metropolis (-Hastings) dengan varian yang lebih kecil dalam distribusi proposal. Dalam keadaan khusus, solusi kedua harus menyatu lebih cepat ke distribusi target. Baik?
AstrOne
Sekarang terkait dengan kode Metropolis-Hastings saya berpikir untuk mengganti R=exp(logR)dengan ini: R=exp(logR)*(dnorm(y[t-1],y.prop,my_sigma)/dnorm(y.prop,y[t-1],my_sigma))untuk MH berjalan acak dan non-acak. Apakah itu benar?
AstrOne
1
Pada dasarnya, tetapi seperti yang saya sebutkan dalam kode metropolis: Anda ingin melakukan perhitungan dalam ruang log. Perhitungan kemungkinan cenderung beroperasi pada nilai yang sangat kecil, sehingga Anda biasanya mendapatkan hasil yang jauh lebih baik dengan menambahkan logaritma dan membuat eksponensial hasil Anda daripada mengalikan nilai mentah secara bersamaan.
David Marx
1
yt-1q(yt|yt-1)=q(yt)yt-1
1
Anda menyatakan "Dalam algoritme berjalan non-acak, distribusi proposal ditetapkan. Dalam varian walk acak, pusat distribusi proposal berubah pada setiap iterasi" yang tidak benar. Versi MH yang tidak berjalan secara acak sebagian besar memiliki proposal yang bergantung pada kondisi rantai Markov saat ini, kadang-kadang bahkan berpusat di kondisi ini. Contoh utama adalah algoritma Langevin MCMC. Ketika proposal diperbaiki, ini adalah algoritma MH independen .
Xi'an