生日问题的模拟计算

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")

unnamed-chunk-1-1.png

模拟的方法

首先构造一个函数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()

unnamed-chunk-14-1.png

记得关闭cluster，由于是通过regiterDoParallel创建的隐式cluster，故通过stopImplicitCluster关闭。

stopImplicitCluster()

可以看到，使用并行计算后，耗时比串行的少。

参考文献

[1]. John Rice. 数理统计与数据分析. 机械工业出版社
[2]. http://michaeljkoontz.weebly.com/uploads/1/9/9/4/19940979/parallel.pdf