生日问题的模拟计算
leengsmile
2016年9月25日
生日问题是概率论中的经典问题，本文介绍R语言模拟生日问题的结果。
概率的方法
生日问题的简单描述是：[1]
若房间中有$n$个人，则至少两个人生日在同一天的概率是多少。
文献[1]中指出，当房间中的人数$n = 23$时，至少有两人的生日在同一天的概率超过0.5。
n <- seq(5, 50, by = 5)
probs <- sapply(n ,
function(t) {
obs <- seq(365, by = -1, length.out = t) / 365
prob <- 1 - prod(obs)
}
)
plot(n, probs, col = "green4")
模拟的方法
首先构造一个函数birthday，计算$num$个人在同一房间中时，至少两人生日相同的概率
require(magrittr)
birthday <- function(num, times = 1000) {
replicate(times,
{
sample(1:days, size = num, replace = TRUE) %>%
table %>%
is_greater_than(1) %>%
any
}
) %>% mean
}
可以模拟当num = 23是对应的概率
days <- 365
birthday(23)
## [1] 0.486
下面计算不同num下的生日问题概率，一个比较简单的版本是使用sapply
system.time(probs_1 <- sapply(n, birthday))
## user system elapsed
## 4.6 0.0 4.8
耗时较长，可以使用foreach包，结合doParallel，做并行计算。
首先引入相应的包
require(foreach)
require(doParallel)
通过registerDoParallel注册foreach的后端，这样才能够做并行的计算。计算结束后，记得停掉相应的后端。
registerDoParallel(cl = 4)
首先来看下面的程序
foreach(i = 1:3) %do% sqrt(i)
## [[1]]
## [1] 1
##
## [[2]]
## [1] 1.414214
##
## [[3]]
## [1] 1.732051
大致流程是foreach通过%do%执行birthday函数。do常用语测试并行语句，如果没有问题，改为dopar就以并行的方式执行语句。
foreach中的num表示循环变量(iteration variable)，可以不止一个循环变量，比如
foreach(a = 1:3, b = 11:14) %do% (a + b)
## [[1]]
## [1] 12
##
## [[2]]
## [1] 14
##
## [[3]]
## [1] 16
foreach的.combine表示已何种方式组合执行的结果，默认的是list，也可以指定为c，表示连接成vector。
foreach(a = 1:3, b = 11:14, .combine = "c") %do% (a + b)
## [1] 12 14 16
也可以指定其他的函数，比如cbind
foreach(i = 1:4, .combine = "cbind") %do% rnorm(5)
## result.1 result.2 result.3 result.4
## [1,] 1.373734919 0.3007412 0.3716988 0.93563782
## [2,] 1.330487498 -0.1006352 -0.6942910 -0.20802592
## [3,] -0.065942480 -0.7285663 0.6705998 -2.40844513
## [4,] 0.009256322 -0.2584815 -1.6157844 -0.01064266
## [5,] -0.293872768 -0.8329582 1.0610386 0.63913900
也可以自定义函数，比如
fun <- function(x, y) 1 / ( 1/x + 1/y)
foreach(i = 1:3, j = 2:4, .combine = fun) %do% (i + j)
## [1] 1.478873
foreach函数知道c、cbind、rbind可以接收很多参数，默认的上限是100。然而对于其他函数，比如+，foreac假定该函数只接受两个参数。如果函数确实接受多个参数，可以将.multicombine设置为TRUE
fun.many <- function(x, ...) c(..., x)
foreach(i = 1:3, .combine = fun.many, .multicombine = TRUE) %do% i
## [1] 2 3 1
执行的结果是1， 2， 3。调用fun.many执行合并时，x对应于1，...对应于2, 3，最后的结果"2, 3"在前，紧跟着"1"。
回到之前的生日问题
system.time(
foreach(num = n, .combine = "c", .packages = c("magrittr")) %dopar%
birthday(num = num) -> ans
)
## user system elapsed
## 0.01 0.00 2.79
.packages是表明dopar后的执行需要用到.packages指定的包里面的函数或者其他。
使用ggplot2绘制相应的模拟
require(ggplot2)
ggplot(data.frame(num = n, prob = ans), aes(x = num, y = prob)) +
geom_line(color = "cyan") + geom_point()
记得关闭cluster，由于是通过regiterDoParallel创建的隐式cluster，故通过stopImplicitCluster关闭。
stopImplicitCluster()
可以看到，使用并行计算后，耗时比串行的少。
参考文献
[1]. John Rice. 数理统计与数据分析. 机械工业出版社
[2]. http://michaeljkoontz.weebly.com/uploads/1/9/9/4/19940979/parallel.pdf
