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

9

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

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

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

q(x)=λ expλx

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

Niech więc p(x) być pdf z XU(0,π), i pozwól Yf(X): celem jest teraz oszacowanie

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

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 powtarzaneN razy, aby uzyskać wiele oszacowań μ, i za każdym razem sprawdza się, czy przedział 95% obejmuje rzeczywistą średnią, czy nie.

Jak widać, dla λ=20faktyczny zasięg wynosi zaledwie 0,19. I rośnieB do wartości takich jak 106nie pomaga (zasięg jest jeszcze mniejszy, 0,15). Dlaczego to się dzieje?

DeltaIV
źródło
1
Użycie nieskończonej funkcji ważności wsparcia dla skończonej całki wsparcia nie jest optymalne, ponieważ część symulacji służy do symulacji zer, że tak powiem. Przynajmniej skróć wykładniczy oπ, co jest łatwe do zrobienia i symulacji.
Xi'an
@ Xi'an, oczywiście, zgadzam się, jeśli musiałbym ocenić tę całkę za pomocą próbkowania ważności, nie użyłbym tego rozkładu ważności, ale próbowałem odpowiedzieć na pierwotne pytanie, które wymagało użycia rozkładu wykładniczego. Mój problem polegał na tym, że nawet jeśli to podejście jest dalekie od optymalnego, zasięg powinien nadal się zwiększać (średnio) asB. I właśnie to pokazał Greenparker.
DeltaIV

Odpowiedzi:

3

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

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

Za pomocą λ=1 Będę w stanie zbadać pełną przestrzeń całkowitą od 0 do πi wydaje się, że tylko kilka losowań πzostaną zmarnowane. Teraz ponownie uruchamiam kod i zmieniam tylkoλ=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 λ, 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=104 do B=106, to tylko przypadkowe zdarzenie, na podstawie tego, którego używasz N=100replikacje. Przedział ufności dla prawdopodobieństwa pokrycia wB=104 jest,

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

Więc nie można tak naprawdę powiedzieć, że rośnie B=106 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=104, prawdopodobieństwo pokrycia wynosi .123 i przy B=106 prawdopodobieństwo pokrycia wynosi .158.

Teraz przedział ufności w okolicach .123 wynosi

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

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

Greenparker
źródło
Tak, wiem, że zasięg zmienia się z λ: w szczególności uzyskano najlepszy zasięg dla 0.1<λ<2. Rozumiem teraz, że ponieważ CI dla średniej próbki opiera się na CLT, jest to wynik asymptotyczny. Dlatego może się tak zmieniaćλwpływa na liczbę próbek potrzebnych do podejścia do „asymptotycznego reżimu”, że tak powiem. Ale chodzi o to, dlaczegoλ=20zasięg maleje od wielkości próby104 do wielkości próbki 106? Z pewnością powinien wzrosnąć, jeśli słaby zasięg był spowodowany tylko wysokimλwartość?
DeltaIV
1
@DeltaIV Dokonałem edycji, aby odpowiedzieć na to pytanie. Istotą jestN=100to za mało replikacji, by powiedzieć coś z całą pewnością.
Greenparker,
1
ah genialne! Nie pomyślałem o utworzeniu przedziału ufności dla samej proporcji pokrycia , a nie tylko dla średniej. Podobnie jak nitpick, nie użyłbym przedziału ufności Walda dla przedziału ufności proporcji. Ponieważ jednak proporcja jest mniejsza od 0 i 1, a liczba powtórzeń wynosi (w drugim przypadkuN=1000) stosunkowo duże, prawdopodobnie użycie interwału Wilsona lub Jeffreysa nie zrobiłoby żadnej różnicy. Poczekam chwilę, aby zobaczyć, czy są inne odpowiedzi, ale powiedziałbym, że w pełni zasługujesz na +100 :)
DeltaIV