编写一个用于硬币抛掷的随机抽样程序。

Question

假设我们有以下情况：

- 有一枚硬币，如果正面朝上，则下一次翻转为正面的概率为0.6（反面朝上则下一次翻转为反面的概率也为0.6） - 一个班级里有100名学生 - 每个学生随机翻转这枚硬币多次 - 第n个学生的最后一次翻转不会影响第n+1个学生的第一次翻转（即当下一个学生翻转硬币时，第一次翻转正反面的概率为0.5，但该学生的下一次翻转取决于前一次翻转结果）

以下是用R代码表示这个问题的方式：

library(tidyverse)

set.seed(123)
ids <- 1:100
student_id <- sort(sample(ids, 100000, replace = TRUE))
coin_result <- character(1000)
coin_result[1] <- sample(c("H", "T"), 1)

for (i in 2:length(coin_result)) {
  if (student_id[i] != student_id[i-1]) {
    coin_result[i] <- sample(c("H", "T"), 1)
  } else if (coin_result[i-1] == "H") {
    coin_result[i] <- sample(c("H", "T"), 1, prob = c(0.6, 0.4))
  } else {
    coin_result[i] <- sample(c("H", "T"), 1, prob = c(0.4, 0.6))
  }
}

my_data <- data.frame(student_id, coin_result)
my_data <- my_data[order(my_data$student_id),]

final <- my_data %>%
    group_by(student_id) %>%
    mutate(flip_number = row_number())
The data looks something like this:

# A tibble: 6 x 3
# Groups:   student_id [1]
  student_id coin_result  flip_number
       <int> <chr>              <int>
1          1 H                      1
2          1 H                      2
3          1 H                      3
4          1 H                      4
5          1 T                      5
6          1 H                      6

我的问题：在这种情况下，假设我对这个硬币没有任何先前的了解（即我只能访问学生提供的数据），我认为这个硬币可能具有“相关概率” - 特别是，我认为上一次抛硬币的结果可能会影响下一次抛硬币的结果。为了测试这个假设，我可以进行以下分析：

• 随机有放回地从学生中抽样，直到你获得与原始数据相同数量的学生。

• 对于每个被选中的学生，随机选择一个起点x和终点y（其中y > x），并选择给定学生在x和y之间的所有可用数据。

• 然后，计算概率和95%置信区间。

• 重复这个过程k次。

以下是我编写上述过程的代码尝试：

library(dplyr)
set.seed(123)

n_boot <- 1000

boot_results2 <- matrix(NA, nrow = n_boot, ncol = 4)
colnames(boot_results2) <- c("P(H|H)", "P(T|H)", "P(H|T)", "P(T|T)")

for (b in 1:n_boot) {

  print(b)
  

  boot_students <- sample(unique(final$student_id), replace = TRUE)
  

  boot_data <- data.frame(student_id = integer(0), coin_result = character(0), stringsAsFactors = FALSE)
  
  for (s in boot_students) {

    student_data <- final %>% filter(student_id == s)
    

    x <- sample(nrow(student_data), 1)
    y <- sample(x:nrow(student_data), 1)
    

    student_data <- student_data[x:y, ]
    

    boot_data <- rbind(boot_data, student_data)
  }
  

  p_hh <- mean(boot_data$coin_result[-1] == "H" & boot_data$coin_result[-nrow(boot_data)] == "H")
  p_th <- mean(boot_data$coin_result[-1] == "H" & boot_data$coin_result[-nrow(boot_data)] == "T")
  p_ht <- mean(boot_data$coin_result[-1] == "T" & boot_data$coin_result[-nrow(boot_data)] == "H")
  p_tt <- mean(boot_data$coin_result[-1] == "T" & boot_data$coin_result[-nrow(boot_data)] == "T")
  
  boot_results2[b, ] <- c(p_hh, p_th, p_ht, p_tt)
}

我的问题：虽然代码似乎在运行，但运行时间很长。我也不确定我是否写对了。

有人可以请指导我如何正确地做吗？

谢谢！

注：可选的代码用于可视化结果：

library(ggplot2)

boot_results_long2 <- as.data.frame(boot_results2)
boot_results_long2$iteration <- 1:n_boot
boot_results_long2 <- boot_results_long2 %>%
  gather(key = "coin", value = "probability", -iteration)


ggplot(boot_results_long2, aes(x = iteration, y = probability, color = coin)) +
  geom_line() +
  labs(x = "Iteration", y = "Probability", color = "Coin") +
  scale_color_discrete(labels = c("P(H|H)", "P(T|H)", "P(H|T)", "P(T|T)"))

- stats_noob

---

嗨 @stats_noob！ - Mark

一些想法：1. 我对统计数据不是百分之百自信，因为这不是我的专业领域——如果你对该观点的统计基础不确定，你可能需要在Cross Validated上确认一下。 - Mark

2. 它运行缓慢的原因很可能是你在使用for循环（而且可能还有嵌套的for循环）。如果你能将代码改写成使用map或apply，那么速度可能会稍微提升一些。你可以在这里或者Code Review上获取更多建议，我不确定。

- Mark

2. 它运行缓慢的原因很可能是你在使用for循环（而且可能还有嵌套的for循环）。如果你能将代码改写成使用map或apply，那么速度可能会有所提升。你可以在这里或者在代码审查网站上获取更多建议，我不确定。

- undefined

Answer 1

看起来主要的瓶颈在于内循环。我们可以通过替换以下内容，将内循环的速度提高约20倍：

tictoc::tic()
for (s in boot_students) {
  student_data <- final %>% filter(student_id == s)
  x <- sample(nrow(student_data), 1)
  y <- sample(x:nrow(student_data), 1)
  student_data <- student_data[x:y, ]
  boot_data <- rbind(boot_data, student_data)
}
tictoc::toc()
# around 2.5s on my machine

使用

tictoc::tic()
boot_data <- final %>%
    left_join(
      final %>%
        ungroup() %>%
        summarize(n = n(), .by = student_id) %>%
        rowwise() %>%
        mutate(x = sample(1:n, 1),
               y = sample(x:n, 1))
    ) %>%
    filter(flip_number >= x, flip_number <= y)
tictoc::toc()
# around 0.1s on my machine

您的原始循环包含了一些低效的步骤：

1. 为每个1000名学生创建单独的最终子集（我们可以跳过这一步）
2. 从x:n中为每个学生随机选择起始点（x）和结束点（y）（让我们只做一次，使用向量化方法）
3. 对数据进行子集操作（让我们只做一次而不是1000次）
4. 将数据附加到之前学生的数据中（如果一开始就不将数据按学生分离，我们可以跳过此步骤）

更高效的做法是先为所有学生执行步骤（2），然后再执行步骤（3），跳过步骤1和4。其中步骤4尤其耗费资源，请参阅R地狱（"Growing Objects"章节）： https://www.burns-stat.com/pages/Tutor/R_inferno.pdf

我相信这样可以进一步加快速度，但也许目前已经足够"快"了。

- Jon Spring

---

非常感谢您的回答！我现在会尝试运行这个代码！最终结果会以相同的格式呈现吗？还能使用ggplot代码制作相同的图表吗？ - stats_noob

我是这么期望的，请你如果发现有什么不对劲的地方，一定要告诉我。 - Jon Spring

@JonSpring：谢谢！我会尝试使用ggplot2来制作一个图表... - stats_noob

Answer 2

如果我们想要对每个student_id进行无权重的抽样，我们可以采用类似@JonSpring的精彩答案的方法，通过分组而不是连接来实现。对于我的计算机来说，使用data.table方法比下面的dplyr方法快大约3倍。

my_sample = function(data) {
  x = sample(nrow(data), 1L)
  y = sample(x:nrow(data), 1L)
  return(data[x:y,])
}

## dplyr
final %>%
  group_by(student_id) %>%
  group_modify(~(my_sample(.x)))

## data.table
library(data.table)
finalDT = as.data.table(ungroup(final))
finalDT[, my_sample(.SD), student_id]

如果我们希望执行类似于 boot_students = sample(unique(final$student_id), replace = TRUE) 的操作，并对它们进行循环，那么使用 data.table 的解决方案应该是相对高效的，因为我们可以事先设置一个键，然后循环遍历所有学生进行筛选。

setkey(finalDT, student_id)
boot_students = sample(unique(final$student_id), replace = TRUE)

lapply(boot_students,
       function (student)  finalDT[student_id == student, my_sample(.SD)]) |>
  rbindlist()

对我来说，它比OP快大约20倍，并且提供与原始内部循环相同的结果。如果我们使用Rcpp，这种方法可能会更高效，因为现在有一些方法可以在Rcpp中对data.table进行子集操作，这应该更有效率。但是...可能已经足够快了；）。

另外，我们可能可以更高效地使用table和potentially prop.table()来完成最后一部分。

## from above
boot_data = lapply(boot_students,
       function (student)  finalDT[student_id == student, my_sample(.SD)]) %>%
  rbindlist()

## table and prop.table on separate lines so intermediate results can be seen
table(head(boot_data$coin_result, -1), tail(boot_data$coin_result, -1)) |>
  prop.table()

##              H         T
##    H 0.2820237 0.2035307
##    T 0.2035307 0.3109150

要绘制图表，我们可以简单地将OP的内部循环替换为上述示例之一。或者，这也可以是使用data.table和ggplot进行绘图的一种方式。与geom_line()不同，它可能不适用于离散数据，这里使用了geom_boxplot()，但可以更深入地考虑如何可视化它。

library(data.table)
library(ggplot2)

## initializing data.table 
finalDT = as.data.table(ungroup(final))
setkey(finalDT, student_id)

## helper methods to subset and then transform into prop.table
my_sample = function(data) {
  x = sample(nrow(data), 1L)
  y = sample(x:nrow(data), 1L)
  return(data[x:y,])
}

create_prop_table = function (x) {
  table(head(x, -1L), tail(x, -1L)) |>
    prop.table()
}

## subset up front for small optimization
unique_students = unique(finalDT[["student_id"]])

## lapply is just a loop
lapply(seq_len(n_boot), 
       function(iteration) {
         boot_students = sample(unique_students, replace = TRUE)
         finalDT[.(boot_students), my_sample(.SD), by = .EACHI
                 ][, as.data.table(create_prop_table(coin_result))
                   ][,c("iteration", "coin") := .(iteration, paste0("P(",V1, "|", V2, ")"))]
       }) |> 
  rbindlist() |> ## combine 100 simulations into one data.table
  ggplot(aes(x = coin, y = N, color = coin)) +
  geom_boxplot() +
  labs(x = "Coin", y = "Probability", color = "Coin") +
  scale_color_discrete(labels = c("P(H|H)", "P(T|H)", "P(H|T)", "P(T|T)")) +
  ggtitle(paste0("Randomized trials based on ", n_boot, " bootstrap simulations"))

- Cole

---

非常感谢您的回答！能否将这些结果可视化呢？ - stats_noob

我的回答末尾的tmp与你的内部循环结束后的boot_data是一样的。可视化结果也是一样的。不过，我已经更新了回答，将tmp重命名为boot_data并提供了一个绘图方法。在我的笔记本电脑上，运行模拟和生成图表大约需要4.5秒，比我的个人电脑上原始内部循环快了大约2秒，而那需要运行100次。 - Cole

@stats_noob 这个回答是否还有什么遗漏的地方？以便被选为最佳答案？ - Cole