Incorrect answer from Monte Carlo integration

Question

How can I approximate the integral of [x^4 * sin(x)]/ [exp(1)^(x/5)] (0 to +Inf) with Monte Carlo method in R?

What I tried to do is

set.seed(666)
func1 <- function(x)
    {
      (x^4 * sin(x))/exp(1)^(x/5)
    }

n <- 1000000

x <- rexp(n, 0.2)
f <- func1(x)
E <- mean(f)

but the result of E is not right

Answer 1

If you're going to sample from exponential, it shouldn't be used again in the function.

From code

set.seed(32345)

func <- function(x) { (x^4 * sin(x)) }

n <- 10000000

x <- rexp(n, 0.2)
f <- func(x)
E <- mean(f)

I'm getting the answer

[1] 13.06643

UPDATE

It fluctuates, and fluctuates badly.

Lest first start with the right answer which according to Mathematica is equal to 4453125/371293 = 11.9936.

I transformed integral from

I = ∫ dx exp(-x/5) x⁴ sin(x)

using substitution y=x/5 to

I = 5⁵ Γ(5) ∫ dy exp(-y) y^5-1 / Γ(5) sin(5*y)

Everything but sin(5*y) is normalized gamma distribution, which we will use to sample, and sin(5*y) will be our function to compute mean value.

And used following trick together with large number of samples: I split calculation of positive values and negative values. It helps if you have fluctuating answer with values canceling each other. I did calculation in batches as well. Gamma function of 5 is just 4! (factorial)

Code

set.seed(32345)

N  <- 10000000 # number of samples per batch
NN <- 640      # number of batches

pos <- rep(0, NN) # positive values
neg <- rep(0, NN) # negative values

for(k in 1:NN) { # loop over batches
    y <- rgamma(N, shape=5, scale=1)
    f <- sin(5.0 * y)
    pnf <- ifelse(f > 0.0, f, 0.0)
    pos[k] <- mean(pnf)
    pnf <- ifelse(f < 0.0, -f, 0.0)
    neg[k] <- mean(pnf)
    print(k)
}

mean(pos)
sd(pos)/sqrt(NN)

mean(neg)
sd(neg)/sqrt(NN)

5*5*5*5*5*4*3*2*(mean(pos) - mean(neg))

Output

> mean(pos)
[1] 0.3183912
> sd(pos)/sqrt(NN)
[1] 4.749269e-06
> 
> mean(neg)
[1] 0.3182223
> sd(neg)/sqrt(NN)
[1] 5.087734e-06
> 
> 5*5*5*5*5*4*3*2*(mean(pos) - mean(neg))
[1] 12.67078

You could see that we really compute difference of two very close values, this is why it is hard to get convergence. It took a bit over 20 minutes to compute on my Xeon workstation.

And with different seed=12345

> mean(pos)
[1] 0.3183917
> sd(pos)/sqrt(NN)
[1] 4.835424e-06
> 
> mean(neg)
[1] 0.3182268
> sd(neg)/sqrt(NN)
[1] 4.633129e-06
> 
> 5*5*5*5*5*4*3*2*(mean(pos) - mean(neg))
[1] 12.36735