Aku mencoba untuk menjawab pertanyaan Evaluasi terpisahkan dengan Pentingnya sampel metode dalam R . Pada dasarnya, pengguna perlu menghitung
menggunakan distribusi eksponensial sebagai distribusi kepentingan
dan temukan nilai yang memberikan perkiraan yang lebih baik untuk integral (itu ). Saya menyusun kembali masalah sebagai evaluasi nilai rata-rata dari atas : integralnya kemudian hanya . self-study
Jadi, misalkan menjadi pdf dari , dan biarkan : tujuannya sekarang adalah memperkirakan
menggunakan sampel penting. Saya melakukan simulasi di R:
# clear the environment and set the seed for reproducibility
rm(list=ls())
gc()
graphics.off()
set.seed(1)
# function to be integrated
f <- function(x){
1 / (cos(x)^2+x^2)
}
# importance sampling
importance.sampling <- function(lambda, f, B){
x <- rexp(B, lambda)
f(x) / dexp(x, lambda)*dunif(x, 0, pi)
}
# mean value of f
mu.num <- integrate(f,0,pi)$value/pi
# initialize code
means <- 0
sigmas <- 0
error <- 0
CI.min <- 0
CI.max <- 0
CI.covers.parameter <- FALSE
# set a value for lambda: we will repeat importance sampling N times to verify
# coverage
N <- 100
lambda <- rep(20,N)
# set the sample size for importance sampling
B <- 10^4
# - estimate the mean value of f using importance sampling, N times
# - compute a confidence interval for the mean each time
# - CI.covers.parameter is set to TRUE if the estimated confidence
# interval contains the mean value computed by integrate, otherwise
# is set to FALSE
j <- 0
for(i in lambda){
I <- importance.sampling(i, f, B)
j <- j + 1
mu <- mean(I)
std <- sd(I)
lower.CB <- mu - 1.96*std/sqrt(B)
upper.CB <- mu + 1.96*std/sqrt(B)
means[j] <- mu
sigmas[j] <- std
error[j] <- abs(mu-mu.num)
CI.min[j] <- lower.CB
CI.max[j] <- upper.CB
CI.covers.parameter[j] <- lower.CB < mu.num & mu.num < upper.CB
}
# build a dataframe in case you want to have a look at the results for each run
df <- data.frame(lambda, means, sigmas, error, CI.min, CI.max, CI.covers.parameter)
# so, what's the coverage?
mean(CI.covers.parameter)
# [1] 0.19
Kode ini pada dasarnya adalah implementasi langsung dari pentingnya pengambilan sampel, mengikuti notasi yang digunakan di sini . Sampling kepentingan kemudian diulang kali untuk mendapatkan beberapa perkiraan , dan setiap kali pemeriksaan dilakukan pada apakah interval 95% mencakup rata-rata aktual atau tidak.
Seperti yang Anda lihat, untuk cakupan sebenarnya hanya 0,19. Dan meningkatkan ke nilai-nilai seperti tidak membantu (cakupannya bahkan lebih kecil, 0,15). Mengapa ini terjadi?
sumber
Jawaban:
Pentingnya pengambilan sampel cukup sensitif terhadap pilihan distribusi yang penting. Karena kamu memilihλ = 20 , sampel yang Anda gunakan 1 / 20 dengan varian 1 / 400 . Ini distribusi yang Anda dapatkan
rexp
akan memiliki rata-rataNamun, integral yang ingin Anda evaluasi berubah dari 0 menjadiπ= 3.14 . Jadi, Anda ingin menggunakan aλ yang memberi Anda kisaran seperti itu. saya menggunakanλ = 1 .
Menggunakanλ = 1 Saya akan dapat menjelajahi ruang integral penuh dari 0 hingga π , dan sepertinya hanya sedikit yang berakhir π akan sia-sia. Sekarang saya jalankan kembali kode Anda, dan hanya berubahλ = 1 .
Jika Anda bermain-main denganλ , Anda akan melihat bahwa jika Anda membuatnya sangat kecil (0,00001) atau besar, probabilitas cakupannya akan buruk.
Sunting -------
Mengenai kemungkinan cakupan berkurang setelah Anda pergiB =104 untuk B =106 , Itu hanya kejadian acak, berdasarkan fakta yang Anda gunakan N= 100 replikasi. Interval kepercayaan untuk probabilitas cakupan diB =104 adalah,
Jadi Anda tidak bisa mengatakan peningkatan ituB =106 secara signifikan menurunkan kemungkinan cakupan.
Bahkan dalam kode Anda untuk seed yang sama, ubahN= 100 untuk N= 1000 , lalu dengan B =104 , probabilitas jangkauan adalah .123 dan dengan B =106 probabilitas cakupan adalah .158 .
Sekarang, interval kepercayaan sekitar 0,123 adalah
Jadi, sekarang denganN= 1000 replikasi, Anda mendapatkan bahwa probabilitas cakupan meningkat secara signifikan.
sumber