Cobertura menor que o esperado para amostragem importante com simulação

9

Eu estava tentando responder à pergunta Avaliar integral com Importância método de amostragem em R . Basicamente, o usuário precisa calcular

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

usando a distribuição exponencial como a distribuição de importância

q(x)=λ expλx

e encontre o valor de que fornece a melhor aproximação à integral (é ). Eu refiz o problema como a avaliação do valor médio de sobre : a integral é então apenas . μ f ( x ) [ 0 , π ] π μλself-studyμf(x)[0,π]πμ

Assim, seja o pdf de e deixe : o objetivo agora é estimarX U ( 0 , π ) Y f ( X )p(x)XU(0,π)Yf(X)

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

usando amostragem de importância. Eu realizei uma simulação em 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

O código é basicamente uma implementação direta de amostragem de importância, seguindo a notação usada aqui . A amostragem de importância é repetida vezes para obter estimativas múltiplas de , e cada vez que é feita uma verificação sobre se o intervalo de 95% cobre a média real ou não.μNμ

Como você pode ver, para a cobertura real é de apenas 0,19. E aumentar para valores como não ajuda (a cobertura é ainda menor, 0,15). Por que isso está acontecendo?B 10 6λ=20B106

DeltaIV
fonte
11
O uso de uma função de importância de suporte infinito para uma integral de suporte finito não é ideal, pois parte das simulações é usada para simular zeros, por assim dizer. Pelo menos trunque o exponencial em , o que é fácil de fazer e simular. π
Xian
@ Xi'an, concordo que, se tivesse que avaliar essa integral por Importance Sampling, não usaria essa distribuição de importância, mas estava tentando responder à pergunta original, que exigia o uso da distribuição exponencial. Meu problema era que, mesmo que essa abordagem esteja longe de ser ótima, a cobertura ainda deve aumentar (em média) como . E foi isso que Greenparker mostrou. B
DeltaIV

Respostas:

3

A amostragem de importância é bastante sensível à escolha da distribuição de importância. Como você escolheu , as amostras que você desenhar usando terão uma média de com variação de . Esta é a distribuição que você recebe1 / 20 1 / 400λ=20rexp1/201/400

insira a descrição da imagem aqui

No entanto, a integral que você deseja avaliar varia de 0 a . Então, você deseja usar um que ofereça esse intervalo. Eu uso .λ λ = 1π=3.14λλ=1

insira a descrição da imagem aqui

Usando , poderei explorar todo o espaço integral de 0 a e parece que apenas alguns desvios sobre serão desperdiçados. Agora, executei novamente seu código e apenas .π π λ = 1λ=1ππλ=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

Se você brincar com , verá que, se você o tornar realmente pequeno (.00001) ou grande, as probabilidades de cobertura serão ruins.λ

EDITAR-------

Em relação à probabilidade de cobertura diminuir depois de passar de para , isso é apenas uma ocorrência aleatória, com base no fato de você usar repetições. O intervalo de confiança para a probabilidade de cobertura em é B = 10 6 N = 100 B = 10 4 0,19 ± 1,96 * B=104B=106N=100B=104

.19±1.96.19(1.19)100=.19±.0769=(.1131,.2669).

Portanto, você não pode realmente dizer que aumentar diminui significativamente a probabilidade de cobertura.B=106

De fato, no seu código para a mesma semente, altere para ; em seguida, com , a probabilidade de cobertura é 0,123 e com probabilidade de cobertura é .N=100N=1000B=104B=106.158

Agora, o intervalo de confiança em torno de .123 é

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

Portanto, agora com repetições, a probabilidade de cobertura aumenta significativamente.N=1000

Greenparker
fonte
Sim, eu sei que a cobertura muda com : em particular, a melhor cobertura é obtida para . Agora, entendo que, como o IC para a média da amostra é baseado no CLT, é um resultado assintótico. Assim, pode ser que a mudança de influencie o número de amostras necessárias para abordar o "regime assintótico", por assim dizer. Mas a questão é: por que com a cobertura diminui do tamanho da amostra para o tamanho da amostra ? Certamente aumentaria, se a cobertura insuficiente fosse apenas devido a um alto valor ? 0,1 < λ < 2 λ λ = 20 10 4 10 6 λλ0,1<λ<2λλ=20104106λ
DeltaIV
11
@ DeltaIV Fiz uma edição para responder a esta pergunta. O ponto principal é que não é replicação suficiente para dizer algo com certeza. N=100
Greenparker
11
ah brilhante! Não pensei em formar o intervalo de confiança para a proporção de cobertura propriamente dita , em vez de apenas para a média. Assim como um nitpick, eu não teria usado o intervalo de confiança de Wald para o intervalo de confiança de uma proporção. No entanto, como a proporção está longe de 0 e 1 e o número de repetições é (no seu segundo caso, ) relativamente grande, provavelmente o uso do intervalo Wilson ou Jeffreys não faria diferença. Vou esperar um pouco para ver se há outras respostas, mas eu diria que você merece plenamente a +100 :)N=1000
DeltaIV