R:随机采样抛硬币组

2024-04-19

我正在使用 R 编程语言。

Suppose:

  • 有一枚硬币,如果它正面朝上,那么下一次抛掷正面的概率是 0.6(如果是反面,那么下一次抛掷反面的概率也是 0.6)
  • 一个班有100名学生
  • 每个学生随机抛掷硬币几次
  • Student_n 的最后一次抛硬币不会影响 Student_n+1 的第一次抛硬币(即当下一个学生抛硬币时,第一次抛硬币正面或反面的概率为 0.5,但该学生的下一次抛硬币取决于上一次抛硬币)

下面是一些代表这个问题的 R 代码:

# https://stackoverflow.com/questions/76192042/r-verifying-the-results-of-coin-flips
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))
  }
}

#tidy up
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())

从这些数据中,我们可以统计 HH、HT、TH 和 TT 序列出现的次数:

head(my_data)
# 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

my_data %>%
  group_by(student_id) %>%
  summarize(Sequence = str_c(coin_result, lead(coin_result)), .groups = 'drop') %>%
  filter(!is.na(Sequence)) %>%
  count(Sequence)

# A tibble: 4 × 2
  Sequence     n
  <chr>    <int>
1 HH       29763
2 HT       19782
3 TH       19775
4 TT       30580

我的问题:使用这些数据,我试图完成以下任务:

  • 步骤 1:使用放回抽样法,选择 100 名学生
  • 步骤 2:对于步骤 1 中选择的给定学生,在其翻转序列中随机选择一个位置(例如,将此位置称为“x”)
  • 步骤 3:对于同一名学生,在翻转序列中随机选择第二个位置(例如,将此位置称为“y”),使得 y > x(即“y”出现在“x”之后)。
  • 步骤 4:对步骤 1 中选择的所有学生重复步骤 2 和步骤 3
  • 步骤5:统计步骤1中选择的所有学生出现HH、HT、TH和TT序列的次数
  • 步骤 6:重复步骤 1 - 步骤 5 多次(例如 1000 次)

例如,假设学生 15 有 6 次翻转:H, H, T, H, T, T- 如果 x = 2 且 y = 5,那么我们将有H, T, H, T

这是我自己解决问题的尝试:

# Set the number of iterations
k <- 1000

# Initialize a data frame 
results <- data.frame(iteration_number = numeric(0),
                      h_given_h = numeric(0),
                      h_given_t = numeric(0),
                      t_given_h = numeric(0),
                      t_given_t = numeric(0))

# Set the number of students to sample
n_students <- length(unique(my_data$student_id))

# Loop over the number of iterations
for (i in 1:k) {
  # Randomly sample student ids with replacement
  sampled_ids <- sample(ids, n_students, replace = TRUE)
  
  # Initialize a data frame to store the sampled data
  sampled_data <- data.frame(student_id = integer(0), coin_result = character(0), stringsAsFactors = FALSE)
  
  # LOOP
  for (j in sampled_ids) {
    # Get data for the current student
    student_data <- my_data[my_data$student_id == j, ]
    
    # Randomly choose a starting and ending point
    x <- sample(nrow(student_data), 1)
    y <- sample(x:nrow(student_data), 1)
    
    # Select the data between the starting and ending point
    selected_data <- student_data[x:y, ]
    
    # Append the selected data to the sampled data frame
    sampled_data <- rbind(sampled_data, selected_data)
  }
  
  final <- sampled_data %>%
    group_by(student_id) %>%
    mutate(flip_number = row_number())
  
  # Calculate the conditional probabilities
  cond_prob <- final %>%
    group_by(student_id) %>%
    summarize(Sequence = str_c(coin_result, lead(coin_result)), .groups = 'drop') %>%
    filter(!is.na(Sequence)) %>%
    count(Sequence) %>%
    mutate(prob = n / sum(n))
  
  # Extract probabilities
  p_HH <- cond_prob$prob[cond_prob$Sequence == "HH"]
  p_HT <- cond_prob$prob[cond_prob$Sequence == "HT"]
  p_TH <- cond_prob$prob[cond_prob$Sequence == "TH"]
  p_TT <- cond_prob$prob[cond_prob$Sequence == "TT"]

  #print(i)
  # Append 
  results[i, ] <- c(i, p_HH, p_HT, p_TH, p_TT)
}

colnames(results) <- c("iteration_number", "h_given_h", "h_given_t", "t_given_h", "t_given_t")

library(ggplot2)

# Convert to long 
results_long <- tidyr::pivot_longer(results, cols = c(h_given_h, h_given_t, t_given_h, t_given_t), names_to = "condition", values_to = "probability")

# Plot
ggplot(results_long, aes(x = iteration_number, y = probability, color = condition)) +
  geom_line() +
  labs(x = "Iteration", y = "Probability", color = "Condition")

我的问题:虽然代码似乎已经运行,但我不确定我是否正确执行了此操作。有人可以帮我确认一下吗?

例如 - HH 和 TT 的行不应该几乎相同......并且 TH 和 HT 的行不应该几乎相同吗?但在我的图表中,情况显然并非如此?在我看来,在给定的迭代中,如果同一个学生在重新采样的数据集中出现 3 次,则第一次的最后一个转换将“泄漏”到第二次的第一次转换中,从而损害结果。

Thanks!


我相信通过重复模拟过程和采样过程,我们可以找到您想要的结果。我将通过创建函数来模拟数据并根据您编写的代码对其进行采样来演示这一点

模拟数据的函数

sim_data <- function(n_students){
ids <- 1:n_students
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))
  }
}

data.frame(student_id, coin_result) %>% 
  arrange(student_id) %>% 
  group_by(student_id) %>%
  mutate(flip_number = row_number())
}

dt <- sim_data(100)

采样功能

Sampling <- function(my_data, n_sample = 100){

  # Step 1
  n_students <- length(unique(my_data$student_id))
  sampled_ids <- sample(1:n_students, n_sample, replace = TRUE)
  sampled_data <- data.frame(student_id = integer(0), 
coin_result = character(0), 
stringsAsFactors = FALSE)
  
  for (j in seq_along(sampled_ids)) {
  # Step 2
    student_data <- my_data[my_data$student_id == sampled_ids[j], ]
    
    x <- sample(nrow(student_data), 1)
    
  # Step 3  
    # There is a edge case where position x is the last flip in the sequence, for 
    # all other cases we shouldn't use position "x" in the sample
    if (x == nrow(student_data)) y <- x else
    y <- sample((x + 1):nrow(student_data), 1)
    selected_data <- student_data[x:y,] %>% 
    mutate(loop_id = j)
    
    sampled_data <- rbind(sampled_data, selected_data)
  }
  
  # Step 5  
  cond_prob <- 
    sampled_data %>%
    group_by(loop_id) %>%
    mutate(flip_number = row_number()) %>% 
    mutate(Sequence = str_c(coin_result, lead(coin_result))) %>%
    filter(!is.na(Sequence)) %>%
    ungroup() %>% 
    count(Sequence) %>%
    mutate(prob = n / sum(n))
  
  
  p_HH <- cond_prob$prob[cond_prob$Sequence == "HH"]
  p_HT <- cond_prob$prob[cond_prob$Sequence == "HT"]
  p_TH <- cond_prob$prob[cond_prob$Sequence == "TH"]
  p_TT <- cond_prob$prob[cond_prob$Sequence == "TT"]
  
  c("h_given_h" = p_HH,"h_given_t" =  p_HT, "t_given_h" = p_TH, "t_given_t" = p_TT)
}

比较各种方法k = 100重复我们得到以下结果:

n_replicate <- 100

results_1 <- 
  replicate(n_replicate,Sampling(dt)) %>% 
  t() %>% 
  as.data.frame() %>%
  rowid_to_column('iteration_number') 

results_2 <- replicate(n_replicate,{
  sim_data(100) %>% 
  Sampling()
  } ) %>% 
  t() %>% 
  as.data.frame() %>%
  rowid_to_column('iteration_number') 

 bind_rows("OP" = results_1,"Sample different data" = results_2, .id = "sample_method") %>% 
   tidyr::pivot_longer(cols = c(h_given_h, h_given_t, t_given_h, t_given_t), names_to = "condition", values_to = "probability") %>%
   summarise(prob = mean(probability), .by = c(sample_method, condition)) %>% 
     pivot_wider(names_from = sample_method, values_from = prob ) 
#> A tibble: 4 × 3
#>  condition    OP `Sample different data`
#>  <chr>     <dbl>                   <dbl>
#> 1 h_given_h 0.283                   0.305
#> 2 h_given_t 0.203                   0.199
#> 3 t_given_h 0.203                   0.199
#> 4 t_given_t 0.311                   0.297

library(ggplot2)

  bind_rows(results_1, results_2, .id = "sample_method") %>% 
    tidyr::pivot_longer(cols = c(h_given_h, h_given_t, t_given_h, t_given_t), names_to = "condition", values_to = "probability") %>% 
    ggplot(aes(x = iteration_number, y = probability, color = condition)) +
  geom_line() +
  labs(x = "Iteration", y = "Probability", color = "Condition") +
    facet_wrap(~sample_method, labeller = labeller(sample_method = c("1" = "OP", "2" = "Sample different data") ) )

Created on 2023-05-11 with reprex v2.0.2 https://reprex.tidyverse.org

Edit

条件概率的计算方式也不理想。在给定“H”的情况下估计“H”概率的最佳方法是计算具有结果“H”的硬币的数量,并在下一个硬币上再次获得结果“H”。您呈现它的方式比 P(H|H) 和 P(T|H) 更接近于估计组合 P(HH) 和 P(HT) 的边际概率。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

R:随机采样抛硬币组 的相关文章

  • R 将多个值与向量进行比较并返回向量[重复]

    这个问题在这里已经有答案了 我有一个向量 A 对于 A 的每个元素 我想检查它是否等于第二个向量 Targets 中的任何元素 我想要一个逻辑值向量 其长度为 A 作为返回 也提到了同样的问题here http r 789695 n4 na
  • 如何在ubuntu的conda环境中更改Rstudio中的R版本

    我在基本系统中安装了 R 4 3 和 Rstudio 在 conda 环境中安装了旧版本的 R 4 2 3 命令which R返回环境中安装的 R 的目录 home 用户 miniconda3 envs anndata2ri pip bin
  • 将第 N 行上的 NA 行插入 data.frames 列表,其中 N 来自列表

    经过几个小时后 我发现自己无法解决以下问题 我有一个数据框列表 我想分别向每个 DF 插入 而不是替换 一行或多行 NA 始终至少一行 要插入的 NA 数量存储在单独的列表中 为了说明这一点 我有以下两个列表 list of datafra
  • 使用字符串中的变量名称访问变量值,R

    Intro 一个数据集有大量的age year变量 age 1990 age 1991 etc 我有一个字符串值数组length age years 表示这些变量 使得age years 1 回报 age 1990 etc Need 我想搜
  • 如何添加链接以从我的 R闪亮应用程序在新窗口中打开 pdf 文件?

    我可以使用 a 从我的 Shiny 应用程序添加到外部站点的超链接 a google href http www google com 但如何创建一个链接来打开 pdf 或类似 文件 看起来应该很简单 但我找不到任何例子 我的问题与此类似
  • 从 n,k 维矩阵数组中减去 n,k 维矩阵

    如果我有一个数组A A lt array 0 c 4 3 5 for i in 1 5 set seed i A i lt matrix rnorm 12 4 3 如果我有矩阵 B set seed 6 B lt matrix rnorm
  • 如何在 R 中合并同名列表中的数据框?

    我有一个包含很多数据框的列表 如果它们具有相同的名称 我想合并它们 即合并所有具有相同名称 a 和 b 的数据框 像这样 a lt aaaaa b lt bbbbb c lt ccccc g lt list df1 lt data fram
  • 为什么数据帧上的 is.vector 不返回 TRUE?

    tl dr R 中的向量到底是什么 长版 R 中很多东西都是向量 例如 数字是长度为 1 的数值向量 is vector 1 1 TRUE 列表也是一个向量 is vector list 1 1 TRUE 好的 所以列表是一个向量 显然 数
  • R中的重叠矩阵

    我有以下数据框 id channel 1 a 1 b 1 c 2 a 2 c 3 a 我想创建并重叠矩阵 它基本上是一个方阵 行和列标签为 a b c 表中的每个条目显示每个通道共有多少个 id 例如 在上面的例子中 矩阵看起来像 a b
  • 列出 R 数据文件的内容而不加载

    我有时用print load myDataFile RData 当我加载数据文件时列出它的内容 有没有办法列出内容而不加载数据文件中包含的对象 我认为如果不加载对象就无法做到这一点 解决方案可能是使用包装器将 R 对象保存到save 该函数
  • 访问或解析 R 中的 summary() 中的元素

    我运行以下 R 命令来进行 Dunnett 测试并获取摘要 如何访问下面线性假设的每一行 这是摘要输出的一部分 基本上我不知道摘要的结构 我尝试使用名称 但它似乎不起作用 因为我没有看到任何命名属性来提供这一点 library multco
  • 如何在 Shiny 中提取动态生成的输入值?

    我正在创建一个闪亮的应用程序 它将根据客户的不同功能为客户生成分数 在我闪亮的应用程序中 我提供了 checkboxGroupInput 来选择所需的功能 根据所选功能 应用程序将动态地将 numericInput 添加到 Web ui 以
  • 多个动态滤镜更新闪亮

    我希望能够让 UI 输入闪亮 并根据用户之前的选择进行自我更新 因此 在下面的示例中 预期的行为是用户选择cyl vsor carb那么这将 过滤数据集mtcars用于创建绘图 即用户根据过滤条件调整绘图并 更新其他过滤器中的剩余输入选择
  • 闪亮的应用程序包:css 和所有 www/ 目录内容

    我正在尝试将 Shiny 应用程序转换为 R 包 但我在处理有关 www 目录以及 松散 文件的所有问题时遇到了问题 我闪亮的应用程序运行得很好 但是当我尝试 打包它 时 它不起作用 我闪亮的应用程序目录 my shiny app R ut
  • 从 data.frame 在 ggplot 图例中添加信息

    我想在图例中添加信息 哪个传感器具有该值 这是我的代码 z lt data frame a c sensor 1 sensor 2 sensor 3 sensor 4 sensor 5 sensor 6 sensor 7 sensor 8
  • R Shinydashboard 自定义 CSS 到 valueBox

    我一直在尝试将 valueBox 的颜色更改为自定义颜色 超出 validColors 中可用的颜色 但一直无法这样做 我知道有一种方法可以使用标签来包含自定义 CSS 但是我无法将它们放在正确的位置 ui lt dashboardPage
  • 在 RMarkdown 输出到 PDF 时缩进而不添加项目符号点或编号

    之前有人问过如何在没有项目符号的情况下缩进文本 RMarkdown 中的点 但这是针对 HTML 输出的 在 RMarkdown 中缩进而不添加项目符号点或数字 https stackoverflow com questions 47087
  • 在ggplot中设置y轴中断

    我在代码中设置中断时遇到困难 我尝试添加breaks seq 0 100 by 20 但似乎无法让它正常工作 本质上我希望 Y 轴从 0 到 100 每 20 个刻度一次 YearlyCI lt read table header T te
  • 增加雷达图中长轴标签的空间

    我想创建一个雷达图ggirahExtra ggRadar 问题是我的标签很长并且被剪掉了 我想我可以通过添加在标签和绘图之间创建更多空间margin margin 0 0 2 0 cm to element text in axis tex
  • 如何为自定义 S3 类实现提取/取子集 ([ [<-, [[ [[<-)] 函数?

    我有一个自定义的 S3 类foo 它在正常的基础上添加了一些自定义行为data frame foo object lt data frame class foo object lt c foo data frame 对于这个类 还应该有一个

随机推荐

  • 将 SVG 从文件加载到画布并取消分组

    我使用 FabricJS 和函数将 SVG 文件上传到画布 fabric loadSVGFromURL url function objects options group fabric util groupSVGElements obje
  • Linux 上共享串口

    我正在使用 Raspberry Pi 进行一个项目 该项目需要能够写入和读取串行端口 但来自不同的程序 程序 A 需要能够从外围设备 A 正在发送数据的串行端口读取数据 程序B需要向串口写入数据 外设B正在监听串口 供参考 本例中程序A是G
  • 将 matlab 中的 find() 转换为 python

    我正在将代码从 Matlab 转换为 Python Matlab中的代码为 x find sEdgepoints gt 0 sNorm lt lowT sEdgepoints x 0 两个数组的大小相同 我基本上是在创建一个掩码 I rea
  • Xcode 14.0 - PackageIndex.findPackages 失败:featureDisabled 警告

    自从我升级到 Xcode 14 0 后 我收到以下警告 PackageIndex findPackages failed featureDisabled 网络搜索没有得到任何结果 我有一个SPM包 但似乎没有任何问题 有人知道如何摆脱这个警
  • 如何在部署的appengine数据库上的eclipse中调试服务器代码?

    我在 Eclipse 中有一个 Google AppEngine Java 项目 我想在 Eclipse 中调试本地代码 但使用 AppEngine 上部署的数据库 到目前为止 我使用带有用户名 密码的远程 API 旧方式 此方法将被弃用
  • 如何获取批处理文件中的字符串长度?

    似乎没有一种简单的方法可以获取批处理文件中字符串的长度 例如 SET MY STRING abcdefg SET A MY STRING LEN 我如何找到字符串的长度MY STRING 如果字符串长度函数处理字符串中所有可能的字符 包括转
  • Chrome 扩展 + 网页视图

    我正在努力寻找这个问题的明确答案 除 Chrome 操作系统外 所有操作系统均已弃用 Chrome 应用 只能在 Chrome 应用中使用 这意味着我不能或不应该在扩展中使用 如果可能 根据进一步的研究 测试和评论 绝对不能在扩展中使用 只
  • postbuild UIAutomation 脚本未在 jenkins 中运行

    我正在尝试做端到端自动化 for an iOS项目 我的目标是自动化持续集成处理与附加UIAutomation脚本作为构建后操作 因此 从用户在 SVN 中检查他的代码开始 直到我们得到自动化测试结果 一切都将是自动化的 Jenkins安装
  • 使用 fb_graph Ruby gem 从 Facebook 检索好友位置

    我正在尝试使用 gem 检索用户所有朋友的位置 fb graph https github com nov fb graph 版本1 7 2 我的权限是 发布流 读取好友列表 离线访问 好友位置 用户位置 我已经对用户进行了身份验证并存储了
  • “不支持”在不指定 RuntimeIdentifier 的情况下构建或发布独立的应用程序

    使用最新的 Visual Studio 2019 我尝试发布 DotNetCore 3 1 WPF 应用程序的 Msix 安装程序 应用程序构建并正确运行 但是当我尝试发布应用程序时出现此错误 It is not supported to
  • 迭代 DFS 与递归 DFS 以及不同的元素顺序

    我编写了一个递归 DFS 算法来遍历图 void Graph
  • eclipse 烦恼:调试和启动工具栏不可用

    我正在运行 Windows XP 和 Eclipse 4 2 2 Build id M20130204 1200 并且我丢失了调试和启动工具栏 我尝试过 Windows gt 重置透视 原始值 和窗口 gt 自定义透视 工具栏可见性和命令组
  • JavaScript 中的错误:对象不是函数

    当我运行下面的代码时 它显示错误object is not a function在控制台中 这个错误就在这一行var todo new Todo contents in my script js文件 我怎样才能让它发挥作用 这是我的 tod
  • 监控网络连接带宽的最佳工具

    我正在寻找一个非常简单的工具来监控所有应用程序的带宽 不需要流量监视等额外功能 我只是对带宽感兴趣 我已经知道 Wireshark 这很棒 但我正在寻找更多类似 TcpView 来自 Sysinternals 的出色工具 以及当前带宽指示的
  • Rails、activerecord 求和然后排序

    我有一个属于用户的工作模型 并且用户有很多工作 我想创建一个 AR 查询来计算每个用户的总工作日数 然后按降序排列 到目前为止 我已经有了这个 但给了我一个错误 列 Job id 必须出现在 GROUP BY 子句中或在聚合函数中使用 wo
  • Phonegap - 在插件委托中从 Objective-c 向 Javascript 发送消息

    我有一个 Phonegap Cordova 插件 在此插件中 我收到来自 javascript 的点击事件 此点击触发使用我的客户端库的文件下载 此文件下载发送事件并调用我的插件中的方法 因为我已将其设置为委托 我无法使用 stringBy
  • java.lang.NoSuchFieldError:没有 Landroidx/compose/foundation/layout/BoxScope$Companion 类型的字段 Companion;

    我是第一次使用 Jetpack Compose 但收到此错误 我还没有弄清楚问题到底出在哪里 但我正在使用单活动架构 如果需要更多信息 请通知我 根据错误信息 问题似乎出在脚手架上 val scaffoldState rememberSca
  • 添加应用程序时 Firebase 数据库被删除

    好的 所以我正在构建一个将在 Play 商店上运行的应用程序 它具有将数据添加到 Firebase 的功能 它无法读取数据 第二个应用程序将保留在我身边 它不会出现在游戏商店中 它用于读取数据 现在我所做的是 假设第一个应用程序有包名称 c
  • 有目的地回到之前的活动

    我有两个活动 当我在第一个活动上按 Enter 时 它将打开第二个活动 它包含一个ListView当我从中选择一个项目时ListView 它将获得其值并返回到第一个活动 这就是我尝试过的 在第二项活动中 listPerasat setOnI
  • R:随机采样抛硬币组

    我正在使用 R 编程语言 Suppose 有一枚硬币 如果它正面朝上 那么下一次抛掷正面的概率是 0 6 如果是反面 那么下一次抛掷反面的概率也是 0 6 一个班有100名学生 每个学生随机抛掷硬币几次 Student n 的最后一次抛硬币