CODEWITHSUNDEEP
rstatisticsintegration
DaveQuinn
DaveQuinn

Reputation: 189

Approximate value for a double integral using monte carlo method in R

I'm trying to implement the monte carlo method in R to find the approximate value for the following double integral:

enter image description here

I have the following code:

mean.estimated <- function(nvals) {
    X <- runif(nvals)
    Y <- runif(nvals)
    sum(exp((2*X + 3*Y)^5))/ nvals
}


monte.carlo <- function(nreps,nvals) {
    estimates <- NULL
    for (i in 1:nreps){
        estimates[i] <- mean.estimated(nvals)
    }
    estimates
}

simvalues <- monte.carlo(200,2000)

But it only produces Inf values. What am i doing wrong?

Upvotes: 0

Views: 3376

Answers (1)

Matthew Lundberg
Matthew Lundberg

Reputation: 42679

With X and Y bounded between 0 and 1, we have that 2*X + 3*Y takes a value between 0 and 5.

When it is the case that 2*X + 3*Y exceeds about 3.72, we have that exp((2*X + 3*Y)^5) is infinite:

> exp((3.72)^5)
[1] Inf

If any one value in the sum is infinite, the sum is infinite. I am not going to compute the odds here, but it is somewhat unlikely that of 2000 samples, every one will have 2*X + 3*Y not exceeding ~3.72. So unlikely that you get Inf for every sample.

Upvotes: 1

Related Questions