r中的生日悖论函数

Question

我是R的初学者，我试图创建一个生日悖论函数，并设法达到了这一点，结果大约是0.5，正如预期的那样。

k <- 23 
sims <- 1000 
event <- 0 
for (i in 1:sims) {
days <- sample(1:365, k, replace = TRUE)
days.unique <- unique(days) 
if (length(days.unique) < k) {
event <- event + 1 } 
answer <- event/sims}
answer

然而，当我试图将其放入函数中时，结果总是0.001。下面是代码:

bdayfunction<- function(k){
sims <- 1000 
event <- 0 
for (i in 1:sims) {
days <- sample(1:365, k, replace = TRUE)
days.unique <- unique(days) 
if (length(days.unique) < k) {
event <- event + 1 } 
answer <- event/sims
return (answer)
}
}

我做错了什么?

Answer 1

您的return不在正确的位置:它在循环中(顺便说一下，您的answer计算也是如此)。

如此:

bdayfunction<- function(k){
sims <- 1000 
event <- 0 
for (i in 1:sims) {
days <- sample(1:365, k, replace = TRUE)
days.unique <- unique(days) 
if (length(days.unique) < k) {
event <- event + 1 }   
}
answer <- event/sims
return (answer)
}

在R中，您可以使用允许您进行分组操作的库。两种主要的是data.table和dplyr。在这里，您可以尝试为所有模拟创建一个长数据帧，然后计算每个模拟的唯一天数，然后计算低于k的发生次数，而不是执行循环。与dplyr:

library(dplyr)
bdayfunction_dplyr <- function(k){  
df <- data.frame(sim = rep(1:sims,each = k),
days = sample(1:365, k*sims, replace = TRUE))
return(
df %>%
group_by(sim) %>%
summarise(plouf = length(unique(days))< k) %>%
summarise(out = sum(plouf)/1000) %>%
pull(out)
)  
}

Indata.table:

library(data.table)
bdayfunction_data.table <- function(k){
dt <- data.table(sim = rep(1:sims,each = k),
days = sample(1:365, k*sims, replace = TRUE))
return(dt[,length(unique(days)),sim][V1<k,.N/1000])
}

您可以测试它们是否提供相同的结果:

set.seed(123)
bdayfunction(23)
[1] 0.515
set.seed(123)
bdayfunction_dplyr(23)
[1] 0.515
set.seed(123)
bdayfunction_data.table(23)
[1] 0.515

现在让我们比较一下速度:

library(microbenchmark)
microbenchmark(initial = bdayfunction(23),
dplyr = bdayfunction_dplyr(23),
data.table = bdayfunction_data.table(23))
Unit: milliseconds
expr     min       lq      mean  median       uq      max neval cld
initial  7.3252  7.56900  8.435564  7.7441  8.15995  24.7681   100  a 
dplyr 12.3488 12.96285 16.846118 13.3777 14.71370 295.6716   100   b
data.table  5.9186  6.24115  6.540183  6.4494  6.75640   8.1466   100  a 

您可以看到，data.table比初始循环稍微快一些，而且写起来更短。