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.