Using mcreplicate Package

In Professor Nahin’s book Will You Be Alive 10 Years From Now?, I read the following interesting problem:

Each of a million men puts his hat into a very large box. Every hat has its owner’s name on it. The box is given a good shaking, and then each men, one after the other randomly draws a hat out of the box. What’s the probability that at least one of the men gets his own hat back?

My immediate response is that I told myself, “I can solve it with Monte Carlo simulation.”

My first version of R code is shown below

set.seed(12345)
N <- 1000

an_exp <- function(n)
{x <- sample(1:n)
 diff <- x - (1:n)
 the_re <- ifelse(0 %in% diff, 1, 0)
 the_re
}

system.time(
  simu_re <- replicate(N, an_exp(n = 1e6))
)
(mean(simu_re))

The code works fine, but it takes nearly three minutes to give me the result. How could I improve on this? After googling, I found that I can use mcreplicate package (https://cran.r-project.org/web/packages/mcreplicate/index.html). So I installed the package and changed my code:

library(mcreplicate)

set.seed(12345)
N <- 1000

an_exp <- function(n)
{x <- sample(1:n)
 diff <- x - (1:n)
 the_re <- ifelse(0 %in% diff, 1, 0)
 the_re
}

system.time(simu_re <- 
              mc_replicate(N, an_exp(1e6), 
                           mc.cores = 2, packages = NULL,
                           varlist = "an_exp", envir = environment())
            )
##    user  system elapsed 
##    0.00    0.00   91.57
(mean(simu_re))
## [1] 0.616

Now the running time is reduced to about 90 seconds. The answer returned by the above R code is 0.62 (if keeping two digits), which agrees with the answer 0.632 in Nahin’s book.

Reference

Nahin, Paul J. (2014). Will You Be Alive 10 Years From Now? p. xvi.

Lingyun Zhang (张凌云)
Lingyun Zhang (张凌云)
Design Analyst

I have research interests in Statistics, applied probability and computation.