Cakupan yang lebih rendah dari yang diharapkan untuk sampel penting dengan simulasi

9

Aku mencoba untuk menjawab pertanyaan Evaluasi terpisahkan dengan Pentingnya sampel metode dalam R . Pada dasarnya, pengguna perlu menghitung

0πf(x)dx=0π1cos(x)2+x2dx

menggunakan distribusi eksponensial sebagai distribusi kepentingan

q(x)=λ exp-λx

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μf(x)[0,π]πμ

Jadi, misalkan menjadi pdf dari , dan biarkan : tujuannya sekarang adalah memperkirakanhal(x)XU(0,π)Yf(X)

μ=E[Y]=E[f(X)]=Rf(x)hal(x)dx=0π1cos(x)2+x21πdx

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.Nμ

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?λ=20B106

DeltaIV
sumber
1
Menggunakan fungsi kepentingan dukungan tak terbatas untuk integral dukungan terbatas tidak optimal karena bagian dari simulasi digunakan untuk mensimulasikan nol, sehingga untuk berbicara. Setidaknya pangkas eksponensial di , yang mudah dilakukan dan disimulasikan. π
Xi'an
@ Xi'an yakin, saya setuju, jika saya harus mengevaluasi integral itu dengan Importance Sampling, saya tidak akan menggunakan distribusi penting itu, tetapi saya mencoba menjawab pertanyaan awal, yang mengharuskan menggunakan distribusi eksponensial. Masalah saya adalah bahwa, bahkan jika pendekatan ini jauh dari optimal, cakupan masih harus meningkat (rata-rata) sebagai . Dan itulah yang ditunjukkan Greenparker. B
DeltaIV

Jawaban:

3

Pentingnya pengambilan sampel cukup sensitif terhadap pilihan distribusi yang penting. Karena kamu memilihλ=20, sampel yang Anda gunakan rexpakan memiliki rata-rata1/20 dengan varian 1/400. Ini distribusi yang Anda dapatkan

masukkan deskripsi gambar di sini

Namun, 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.

masukkan deskripsi gambar di sini

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.

# 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(1,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] .95

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 pergi B=104 untuk B=106, Itu hanya kejadian acak, berdasarkan fakta yang Anda gunakan N=100replikasi. Interval kepercayaan untuk probabilitas cakupan diB=104 adalah,

.19±1.96.19(1-.19)100=.19±0,069=(.1131,.2669).

Jadi Anda tidak bisa mengatakan peningkatan itu B=106 secara signifikan menurunkan kemungkinan cakupan.

Bahkan dalam kode Anda untuk seed yang sama, ubah N=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

.123±1.96.123(1-.123)1000=.123±0,0203=(.102,.143).

Jadi, sekarang dengan N=1000 replikasi, Anda mendapatkan bahwa probabilitas cakupan meningkat secara signifikan.

Greenparker
sumber
Ya, saya tahu bahwa cakupannya berubah dengan λ: khususnya cakupan terbaik diperoleh untuk 0,1<λ<2. Sekarang, saya mengerti bahwa karena CI untuk mean sampel didasarkan pada CLT, ini merupakan hasil asimptotik. Jadi, mungkin saja perubahan ituλmempengaruhi jumlah sampel yang diperlukan untuk mendekati "rezim asimptotik", bisa dikatakan. Tapi intinya, mengapa denganλ=20cakupan berkurang dari ukuran sampel104 untuk ukuran sampel 106? Tentunya itu harus meningkat, jika cakupannya buruk hanya karena tingginyaλnilai?
DeltaIV
1
@DeltaIV Saya mengedit untuk menjawab pertanyaan ini. Intinya adalah,N=100tidak cukup replikasi untuk mengatakan sesuatu dengan pasti.
Greenparker
1
ah brilian! Saya tidak berpikir untuk membentuk interval kepercayaan untuk proporsi cakupan itu sendiri , bukan hanya untuk rata-rata. Sama seperti nitpick, saya tidak akan menggunakan interval kepercayaan Wald untuk interval kepercayaan proporsi. Namun, karena proporsinya jauh dari 0 dan 1 dan jumlah ulangannya (dalam kasus kedua Anda,N=1000) relatif besar, mungkin menggunakan interval Wilson atau Jeffreys tidak akan membuat perbedaan. Saya akan menunggu sedikit untuk melihat apakah ada jawaban lain, tapi saya akan mengatakan Anda sepenuhnya layak +100 :)
DeltaIV