Monte Carlo Integration in R from normal distribution

Question

Monte Carlo Question in R. I have the following equation.

h(x)=x/((3^x)-1)

I am trying to compute a estimator for the integral where random variables are generated from N(0,1).

The following code works in generating the correct value:

set.seed(1)
x<-runif(1000,0,1)
x2<-(x/((3^x)-1))

result<-mean(x2)
print(result)

However this is using uniform r.v's. When I modify to runif to rnorm I am not getting the desired value. I tried multiplying the h(x) equation by the pdf of norm as well.

Attempt at normal:

    set.seed(1)
xlist<-rnorm(n=10000, mean=0, sd=1)
result<-0
for (x in xlist)
{
  a1<-(1*sqrt(2*pi))/((exp((-1/2)*((x-0)/1)^2)))
  a2<-(x/((3^x)-1))
  a3<-a1*a2
  result<-result+a3
}
print (result/n)

Answer 1

Here's how I would go about this. First, let's examine the function we want to integrate:

f <- function(x) (x/((3^x)-1))

curve(f)

We can see that in the interval 0, 1, the maximum and minimum values of the function are between 0 and 1. We can see that the integral is less than 1, and must be some proportion of the unit square bounded by x = 0, x = 1, y = 0 and y = 1.

This means that if we create uniform random points on the unit square, the proportion that fall below the line will approximate the integral of the function. We can simulate the points like this:

set.seed(1)
xvals <- runif(1e6)
yvals <- runif(1e6)

Note that if, for some reason, you are constrained to use rnorm rather than runif, you can get the same effect by doing

set.seed(1)
xvals <- pnorm(rnorm(1e6))
yvals <- pnorm(rnorm(1e6))

In either case our answer is just the proportion of these points that are below our line:

sum(yvals < f(xvals)) / length(xvals)
#> [1] 0.691146

This is pretty close to the actual value of 0.690395, as you can see on Wolfram Alpha.

^{Created on 2022-03-13 by the reprex package (v2.0.1)}