Niższy niż oczekiwano zasięg ważnego próbkowania z symulacją

Starałem się odpowiedzieć na pytanie Ocenić integralny z Znaczenie metody pobierania próbek na badania . Zasadniczo użytkownik musi obliczyć

wykorzystanie rozkładu wykładniczego jako rozkładu ważności

i znajdź wartość $\lambda$co daje lepsze przybliżenie całki (jej self-study). Przekształcam problem jako ocenę wartości średniej$\mu$ z $f(x)$ nad $[0,\pi]$: całka jest wtedy sprawiedliwa $\pi\mu$.

Niech więc $p(x)$ być pdf z $X\sim\mathcal{U}(0,\pi)$, i pozwól $Y\sim f(X)$: celem jest teraz oszacowanie

przy użyciu próbkowania ważności. Przeprowadziłem symulację w języku 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

Kod jest w zasadzie prostą implementacją ważności próbkowania, zgodnie z zastosowaną tutaj notacją . Próbkowanie według ważności jest następnie powtarzane$N$ razy, aby uzyskać wiele oszacowań $\mu$, i za każdym razem sprawdza się, czy przedział 95% obejmuje rzeczywistą średnią, czy nie.

Jak widać, dla $\lambda=20$faktyczny zasięg wynosi zaledwie 0,19. I rośnie$B$ do wartości takich jak $10^6$nie pomaga (zasięg jest jeszcze mniejszy, 0,15). Dlaczego to się dzieje?

Ważność próbkowania jest dość wrażliwa na wybór rozkładu ważności. Ponieważ wybrałeś$\lambda = 20$, próbki, które narysujesz, rexpbędą miały średnią$1/20$ z wariancją $1/400$. To jest dystrybucja, którą otrzymujesz

Jednak całka, którą chcesz ocenić, zmienia się od 0 do $\pi =3.14$. Więc chcesz użyć$\lambda$co daje ci taki zasięg. używam$\lambda = 1$.

Za pomocą $\lambda = 1$ Będę w stanie zbadać pełną przestrzeń całkowitą od 0 do $\pi$i wydaje się, że tylko kilka losowań $\pi$zostaną zmarnowane. Teraz ponownie uruchamiam kod i zmieniam tylko$\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

Jeśli bawisz się z $\lambda$, zobaczysz, że jeśli zrobisz to naprawdę małe (.00001) lub duże, prawdopodobieństwo pokrycia będzie złe.

EDYTOWAĆ-------

Jeśli chodzi o prawdopodobieństwo pokrycia, maleje ono po przejściu $B = 10^4$ do $B = 10^6$, to tylko przypadkowe zdarzenie, na podstawie tego, którego używasz $N = 100$replikacje. Przedział ufności dla prawdopodobieństwa pokrycia w$B = 10^4$ jest,

Więc nie można tak naprawdę powiedzieć, że rośnie $B = 10^6$ znacznie obniża prawdopodobieństwo pokrycia.

W rzeczywistości w swoim kodzie dla tego samego materiału siewnego zmień $N = 100$ do $N = 1000$, a następnie z $B = 10^4$, prawdopodobieństwo pokrycia wynosi .123 i przy $B = 10^6$ prawdopodobieństwo pokrycia wynosi $.158$.

Teraz przedział ufności w okolicach .123 wynosi

Tak więc teraz z $N = 1000$ replikacje, otrzymujesz, że prawdopodobieństwo pokrycia znacznie wzrasta.