0

我想创建一个具有三个参数的函数：样本大小、样本数量、伯努利试验中成功的真实 p。

该函数将输出以下结果：p 的平均估计值（即每个样本的 p 估计值的平均值）、真实标准差的平均估计值（即每个样本的 sd 估计值的平均值），最后是分数95% CI 包含真实 p 的样本的数量。

我想出了以下代码：

   draw_function <- function(si=1, proba= 0.5){
   sample_draw <- rbinom(n= si, 1, prob= proba)
   mean_estimate <-  mean(sample_draw)
   s_estimate <- sqrt(mean_estimate*(1-mean_estimate))
  confidence_interval <- list(c(mean_estimate - 
  1.96*s_estimate/sqrt(si), mean_estimate + 1.96*s_estimate/sqrt(si)))
    list(mean_sample = mean_estimate, s_sample = s_estimate, c_sample = 
  confidence_interval)
   }



   draw_sample <- function(s=1, r=1, p=0.5) {
   for (i in 1:r) {
     inter <-  list(mean_p=c(), mean_s = c(), mean_c = c() )
     samp <- draw_function(si=s, proba= p)
    inter$mean_p = c(inter$mean_p, samp$mean_sample)
  inter$mean_s = c(inter$mean_s, samp$s_sample)
  inter$mean_c <- c(inter$mean_c, samp$c_sample)
    if ( inter[[3]][1] > p & inter[[3]][2] < p ) {
   mean_c <- (mean_c + 1)/i
   } 
  else {
  meanc <- (mean_c + 0)/i
  }
    }

  return(list(mean(inter$mean_p), mean(inter$mean_s), inter$mean_c))

但是，即使在修改了引起我注意的错误之后，此代码也不起作用。

我不断收到此错误：

draw_sample 中的错误（s = 30，r = 1000，p = 0.05）：（
列表）对象不能被强制输入“双”

因此，我很想寻求您的帮助以找到问题并构建这样的功能！谢谢！

4

1 回答 1

0

错误发生在以下部分：

inter[[3]][1] > p & inter[[3]][2] < p

使用 browser() 或逐行运行代码，您会注意到

confidence_interval <- list(c(mean_estimate - 1.96*s_estimate/sqrt(si), 
  mean_estimate + 1.96*s_estimate/sqrt(si)))

返回一个列表本身。因此，如果您想自己获取零件，则必须分别是inter[[3]][[i]][1]和inter[[3]][[i]][2]。

可以进行一些优化，我最近可能会编辑此评论并提出一些建议。:-)

:::Edit::: 一个更深入的答案。R 中的r[dens]函数有一些与之相关的开销。作为加快代码速度的一种简单方法，就是一次运行所有模拟，并以巧妙的方式对子集执行计算。一种方法（如果模拟在内存限制内）如下所示，将所有内容插入矩阵，然后使用 rowMeans 快速计算必要的信息。

draw_sample <- function(s = 1, r = 1, p = .5){
  samples <- matrix(rbinom(n = r * s, size = 1, p = p), ncol = s, nrow = r, byrow = TRUE)
  mean_p <- rowMeans(samples)
  mean_s <- mean_p * (1 - mean_p)
  bound = 1.96 * mean_s / sqrt(s)
  mean_c <- mean(ifelse(mean_p - bound < p & mean_p + bound > p, 1, 0))
  list(mean_p = mean(mean_p), mean_s = mean(mean_s), mean_c = mean_c)
}
于 2019-01-28T18:12:43.433 回答