Eu estava tentando responder à pergunta Avaliar integral com Importância método de amostragem em R . Basicamente, o usuário precisa calcular
usando a distribuição exponencial como a distribuição de importância
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
Assim, seja o pdf de e deixe : o objetivo agora é estimarX ∼ U ( 0 , π ) Y ∼ f ( X )
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.μ
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
fonte
Respostas:
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λ = 20 1 1 / 20 1 1 / 400
rexp
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
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
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 = 104 B = 106 N= 100 B= 104
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= 100 N= 1000 B = 104 B = 106 .158
Agora, o intervalo de confiança em torno de .123 é
Portanto, agora com repetições, a probabilidade de cobertura aumenta significativamente.N= 1000
fonte