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.