假设我们有以下情况:
- 有一个硬币,如果它是正面的,那么下一个硬币是正面的概率是0.6(如果是反面的,那么下一个硬币是反面的概率也是0.6)
- 一个班有100个学生
- 每个学生随机抛硬币几次
- student_n的最后翻转不影响student_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)"))
型
2条答案
按热度按时间bnl4lu3b1#
看起来主要的瓶颈是内部循环。我们可以通过替换以下内容使内部循环的速度提高20倍:
字符串
与
型
你的原始循环包括一些低效的步骤:
1.为1000名学生中的每一个创建单独的final子集(我们可以跳过此操作)
1.为每个学生从x:n中随机选取开始(x)和结束(y)。(让我们做一次,成为vectorized)
1.对数据进行子集化(让我们只做一次,而不是1000x)
1.附加到以前的学生数据(如果我们从不按学生分开数据,我们可以跳过这一步)
更有效率的做法是一次对所有学生执行(2),然后跳过(3),跳过1+4。4是特别昂贵的,见第2章(“成长的对象”)的R地狱:https://www.burns-stat.com/pages/Tutor/R_inferno.pdf
我相信这可以进一步加快,但也许这进入了“目前足够快”的区域。
2izufjch2#
如果我们想在没有任何权重的情况下对每个
student_id
进行采样,我们可以通过分组而不是加入来接近它,类似于@JonSpring的精彩答案。对于我的计算机,data.table
方法比下面的dplyr
方法快3倍。字符串
如果我们想做一些类似于
boot_students = sample(unique(final$student_id), replace = TRUE)
的事情并循环遍历它们,那么data.table
解决方案应该相对性能更好,因为我们可以预先设置一个键,然后循环遍历所有学生以进行过滤。型
对我来说,它比OP快20倍,并提供与原始内部循环相同的结果。如果我们放弃使用
Rcpp
,这种方法可以变得更高效,因为现在有一些方法可以在Rcpp
中子集data.table
,这应该更有效。但是...可能已经够快了).顺便说一句,我们可能可以使用
table
和潜在的prop.table()
更有效地完成最后部分。型
为了绘制图表,我们可以简单地用上面的例子之一替换OP的内部循环。或者,这可以是使用data.table和ggplot进行绘图的方法。这里不是
geom_line()
,它可能不是离散数据的最佳选择,而是geom_boxplot()
,但可以更多地考虑如何将其可视化。型
x1c 0d1x的数据