Geringere Abdeckung als erwartet für die Stichprobenerhebung mit Simulation

Ich habe versucht, die Frage Integral mit Wichtigkeit Stichprobenmethode in R bewerten zu beantworten . Grundsätzlich muss der Benutzer berechnen

Verwenden der Exponentialverteilung als Wichtigkeitsverteilung

und finde den Wert von der die bessere Annäherung an das Integral (es ist ) ergibt . Ich formuliere das Problem als die Bewertung des Mittelwerts von über : Das Integral ist dann nur . $\lambda$self-study$\mu$$f(x)$$[0,\pi]$$\pi\mu$

Sei also das PDF von und sei : Das Ziel ist nun die Schätzung$p(x)$$X\sim\mathcal{U}(0,\pi)$$Y\sim f(X)$

unter Verwendung von Wichtigkeitsstichproben. Ich habe eine Simulation in R durchgeführt:

# 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

Der Code ist im Grunde eine einfache Implementierung der Wichtigkeitsabtastung gemäß der hier verwendeten Notation . Die Wichtigkeitsabtastung wird dann mal wiederholt , um mehrere Schätzungen von , und jedes Mal wird überprüft, ob das 95% -Intervall den tatsächlichen Mittelwert abdeckt oder nicht.$N$$\mu$

Wie Sie sehen können, beträgt die tatsächliche Abdeckung für nur 0,19. Und das Erhöhen von auf Werte wie hilft nicht (die Abdeckung ist mit 0,15 sogar noch kleiner). Warum passiert dies?$\lambda=20$$B$$10^6$

Die Stichprobenerhebung ist sehr empfindlich gegenüber der Wahl der Wichtigkeitsverteilung. Da Sie , haben die Stichproben, mit denen Sie zeichnen, einen Mittelwert von mit einer Varianz von . Dies ist die Verteilung, die Sie erhalten$\lambda = 20$rexp$1/20$$1/400$

Das Integral, das Sie auswerten möchten, reicht jedoch von 0 bis . Sie möchten also ein , das Ihnen einen solchen Bereich bietet. Ich benutze .$\pi =3.14$$\lambda$$\lambda = 1$

Mit ich den gesamten Integralraum von 0 bis erkunden und es scheint, als würden nur ein paar Draws über verschwendet. Jetzt führe ich Ihren Code erneut aus und ändere nur .$\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

Wenn Sie mit , werden Sie feststellen, dass die Abdeckungswahrscheinlichkeiten schlecht sind, wenn Sie es wirklich klein (.00001) oder groß machen.$\lambda$

BEARBEITEN-------

Wenn die Abdeckungswahrscheinlichkeit abnimmt, wenn Sie von auf , ist dies nur ein zufälliges Ereignis, basierend auf der Tatsache, dass Sie Replikationen verwenden. Das Konfidenzintervall für die Abdeckungswahrscheinlichkeit bei beträgt $B = 10^4$$B = 10^6$$N = 100$$B = 10^4$

Man kann also nicht wirklich sagen, dass eine Erhöhung von die Abdeckungswahrscheinlichkeit signifikant senkt.$B = 10^6$

In der Tat im Code für die gleiche Saat, ändert bis , dann mit ist Deckungswahrscheinlichkeit .123 und mit Deckungswahrscheinlichkeit ist .N = 1000 B = 10 4 B = 10 6$N = 100$$N = 1000$$B = 10^4$$B = 10^6$$.158$

Das Konfidenzintervall um .123 beträgt nun

Somit erhalten Sie jetzt mit Replikationen, dass die Abdeckungswahrscheinlichkeit signifikant zunimmt.$N = 1000$