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

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

utiliser la distribution exponentielle comme distribution d'importance

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 , π ] π $\lambda$self-study$\mu$$f(x)$$[0,\pi]$ : l'intégrale est alors juste . $\pi\mu$

Ainsi, soit le pdf deX ∼ U ( 0 , π ) Y ∼ f ( X $p(x)$$X\sim\mathcal{U}(0,\pi)$ , et soit : le but est maintenant d'estimer$Y\sim f(X)$

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.$\mu$

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 $\lambda=20$$B$$10^6$

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 /$\lambda = 20$rexp$1/20$$1/400$

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

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 .$\lambda = 1$$\pi$$\pi$$\lambda = 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$\lambda$ , 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 = 10^4$$B = 10^6$$N = 100$$B = 10^4$

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

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$N = 100$$N = 1000$$B = 10^4$$B = 10^6$$.158$

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

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