R: monte carlo integration using Importance Sampling

pdf-generation, probability, r, statistics

Solution

Generally your approach seems to be correct, but you have to be more careful with the domain over which you want to integrate. In your original example, about 20% of values `rexp(1000, 1.5)` are above 1. The function `dexp(x, rate=1.5)` is not a density function on the interval [0,1]. You have to divide by `pexp(1, rate=1.5)`. So here is what I would do for the importance sampling example:

#Importance sampling Monte Carlo
w <- function(x) dunif(x, 0.01, 1)/dexp(x,rate=1.5) * pexp(1, rate=1.5)
f <- function(x) x^(-0.5)
X <- rexp(1000,rate=1.5)
X <- X[X<=1]
Y <- w(X)*f(X)
c(mean(Y), var(Y))

In your second example the same thing causes the problem. You get negative X and therefore get NA values for log(X). Furthermore, your normal function should be centered at 0.5 with less variance. Here's my approach:

#Without Importance Sampling
set.seed(1909)
X <- runif(1000,0.01,1)
Y <- (1+sinh(2*X)*log(X))^(-1)
c(mean(Y), var(Y))

#Importance sampling Monte Carlo
w <- function(x) dunif(x, 0.01, 1)/dnorm(x, mean=0.5, sd=0.25) * (1-2*pnorm(0, mean=0.5, sd=0.25))
f <- function(x) (1+sinh(2*x)*log(x))^(-1)
X <- rnorm(1000, mean=0.5, sd=0.25)
Y1 <- w(X)
Y2 <- f(X)
Y <- Y1*Y2
Y <- Y[!(is.na(Y2)&Y1==0)]
c(mean(Y), var(Y))

In your second example, I don't really understand what `y` is. Is it just a constant? Then perhaps a Weibull distribution may work.

EDIT: Regarding your additional questions in the comments. (1) Any probability density function should integrate to 1. Therefore `dexp(x, rate=1.5)` is not a density function on the interval [0,1], it only integrates to `pexp(1, rate=1.5)`. However, the function

dexp01 <- function(x, rate){
  dexp(x, rate=rate)/pexp(1, rate=rate)
}

actually integrates to 1:

integrate(dexp, 0, 1, rate=1.5)
integrate(dexp01, 0, 1, rate=1.5)

That's the rationale of including the probability distribution function. If you have a different interval, e.g. [0.3,8], you have to adjust the function accordingly:

dexp0.3_8 <- function(x, rate){
  dexp(x, rate=rate)/(pexp(8, rate=rate)-pexp(0.3, rate=rate))
}
integrate(dexp0.3_8, 0.3, 8, rate=1.5)

(2) Here I choose the variance so that approximately 95% of the values in `rnorm(1000, .5, .25)` were in the interval from 0 to 1 (having many values outside this interval would certainly increase the variance). However, I am not certain that this is the best choice of distribution function. The selection of the importance function is a problem that I am not very familiar with. You could ask on CrossValidated. Same goes for your next question.

Problem

Original source