Couverture plus faible que prévu pour l'échantillonnage d'importance avec simulation

9

Je suis en train de répondre à la question Évaluer une seule pièce avec méthode d' échantillonnage importance dans R . Fondamentalement, l'utilisateur doit calculer

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

utiliser la distribution exponentielle comme distribution d'importance

q(x)=λ expλx

et trouver la valeur de qui donne la meilleure approximation de l'intégrale (c'est ). Je refonte le problème en évaluant la valeur moyenne de surμ f ( x ) [ 0 , π ] π μλself-studyμf(x)[0,π] : l'intégrale est alors juste . πμ

Ainsi, soit le pdf deX U ( 0 , π ) Y f ( X )p(x)XU(0,π) , et soit : le but est maintenant d'estimerYf(X)

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

en utilisant l'échantillonnage d'importance. J'ai effectué une simulation en 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

Le code est fondamentalement une implémentation simple de l'échantillonnage d'importance, suivant la notation utilisée ici . L'échantillonnage d'importance est ensuite répétéμN fois pour obtenir plusieurs estimations de , et chaque fois qu'un contrôle est effectué pour savoir si l'intervalle de 95% couvre la moyenne réelle ou non.μ

Comme vous pouvez le voir, pour la couverture réelle n'est que de 0,19. Et l'augmentation de à des valeurs telles que n'aide pas (la couverture est encore plus petite, 0,15). Pourquoi cela arrive-t-il?B 10 6λ=20B106

DeltaIV
la source
1
L'utilisation d'une fonction d'importance de support infinie pour une intégrale de support fini n'est pas optimale car une partie des simulations est utilisée pour simuler des zéros, pour ainsi dire. Tronquez au moins l'exponentielle en , ce qui est facile à faire et à simuler. π
Xi'an
@ Xi'an bien sûr, je suis d'accord, si je devais évaluer cette intégrale par Importance Sampling, je n'utiliserais pas cette distribution d'importance, mais j'essayais de répondre à la question d'origine, qui nécessitait d'utiliser la distribution exponentielle. Mon problème était que, même si cette approche est loin d'être optimale, la couverture devrait quand même augmenter (en moyenne) de . Et c'est ce que Greenparker a montré. B
DeltaIV

Réponses:

3

L'échantillonnage d'importance est assez sensible au choix de la distribution d'importance. Puisque vous avez choisi , les échantillons que vous dessinez en utilisant auront une moyenne de avec la variance . Voici la distribution que vous obtenez1 / 20 1 / 400λ=20rexp1/201/400

entrez la description de l'image ici

Cependant, l'intégrale que vous souhaitez évaluer passe de 0 à . Vous voulez donc utiliser un qui vous donne une telle plage. J'utilise .π=3.14λλ=1

entrez la description de l'image ici

En utilisant je pourrai explorer l'espace intégral complet de 0 à , et il semble que seuls quelques tirages sur seront gaspillés. Maintenant, je réexécute votre code et je ne change que .λ=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

Si vous jouez avecλ , vous verrez que si vous le rendez vraiment petit (.00001) ou grand, les probabilités de couverture seront mauvaises.

ÉDITER-------

En ce qui concerne la probabilité de couverture diminuant une fois que vous passez de à , il s'agit simplement d'une occurrence aléatoire, basée sur le fait que vous utilisez réplications. L'intervalle de confiance pour la probabilité de couverture à est de B=104B=106N=100B=104

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

Vous ne pouvez donc pas vraiment dire que l'augmentation de diminue considérablement la probabilité de couverture.B=106

En fait, dans votre code pour la même graine, changez en , puis avec , la probabilité de couverture est de 0,123 et avec la probabilité de couverture est de .N = 1000 B = 10 4 B = 10 6 .158N=100N=1000B=104B=106.158

Maintenant, l'intervalle de confiance autour de .123 est

.123±1.96.123(1.123)1000=.123±.0203=(.102,.143).

Ainsi, maintenant avec réplications, vous obtenez que la probabilité de couverture augmente considérablement.N=1000

Greenparker
la source
Oui, je sais que la couverture change avec : en particulier la meilleure couverture est obtenue pour . Maintenant, je comprends que puisque l'IC de la moyenne de l'échantillon est basé sur le CLT, c'est un résultat asymptotique. Ainsi, il se pourrait bien que le changement de influence le nombre d'échantillons nécessaires pour approcher le "régime asymptotique", pour ainsi dire. Mais le fait est que, avec la couverture diminue de la taille de l'échantillon à la taille de l'échantillon ? Elle devrait sûrement augmenter, si une mauvaise couverture n'était due qu'à une valeur élevée? 0,1 < λ < 2 λ λ = 20 10 4 10 6 λλ0.1<λ<2λλ=20104106λ
DeltaIV
1
@DeltaIV J'ai fait une modification pour répondre à cette question. L'essentiel est que n'est pas assez de répétitions pour dire quoi que ce soit avec certitude. N=100
Greenparker
1
ah génial! Je n'ai pas pensé à former l'intervalle de confiance pour la proportion de couverture elle - même , plutôt que juste pour la moyenne. Tout comme une piqûre, je n'aurais pas utilisé l'intervalle de confiance de Wald pour l'intervalle de confiance d'une proportion. Cependant, comme la proportion est éloignée de 0 et 1 et que le nombre de répétitions est (dans votre deuxième cas, ) relativement important, l'utilisation de l'intervalle Wilson ou Jeffreys n'aurait probablement fait aucune différence. J'attendrai un peu pour voir s'il y a d'autres réponses, mais je dirais que vous méritez pleinement le +100 :)N=1000
DeltaIV