通过从 R 中每个试验的概率计算累积概率来学习使用循环

Question

我已经独自挣扎了很长时间才能找到答案。我保证我会尝试从解决方案中学习。为了学习，我想了解如何使用显式循环来做到这一点，但如果您想分享矢量化方法作为奖励，也非常感谢。

假设我每天要玩一次游戏，并且我知道每天获胜的概率。我想要一个函数，它采用概率向量并返回至少一天成功的累积概率。所以如果我连续玩了 3 天并且每天获胜的概率是 0.5，那么我的函数应该返回 "0.875, 0.75, 0.5"

这是我最近一次编写此函数的失败尝试：

prob_cum <- function(prob_today) {
  p_cum <- rep(0, length(prob_today))
  for (i in 1:length(prob_today)) {
    for (j in i:length(prob_today)) {
      p_cum[j] <- p_cum[j-1] - ((1 - p_cum[j-1]) * prob_today[j])
    }
  }
  p_cum
}

prob_daily <- c(.5,.5,.5)
prob_cum(prob_daily)

Answer 1

>  1 - cumprod( 1- c(0.5,0.5,0.5) )
[1] 0.500 0.750 0.875
 # (1- prob_success) is the prob_non_success vector

如果需要，可以轻松包装到函数中。您的初始测试不是一个好的测试，因为它没有透露我在 cumprod 参数中没有从 1 中减去成功向量的原始错误。

 vec<-runif(100)
 prob_cum <- function(prob_today) {
   p_cum <- rep(0, length(prob_today))
   p_cum[1] <- prob_today[1]
   for (j in seq_along(prob_today)[-1]) {
     p_cum[j] <- p_cum[j-1] + ((1 - p_cum[j-1]) * prob_today[j])
   }
   p_cum
 }
 Prob_vec <- function(vec) 1 - cumprod( 1- vec) 
 require(rbenchmark)
 benchmark( prob_cum(vec) , Prob_vec(vec) ,replications=1000)
#           test replications elapsed relative user.self sys.self user.child sys.child
#1 prob_cum(vec)         1000   0.538   59.778     0.532    0.008          0         0
#2 Prob_vec(vec)         1000   0.009    1.000     0.008    0.002          0         0

Answer 2

一次解决每个问题：

你有一个循环 i 没有做任何事情；它只是多次执行相同的计算，并且每次都会覆盖结果（具有相同的结果）。放下那个。

prob_cum <- function(prob_today) {
  p_cum <- rep(0, length(prob_today))
  for (j in i:length(prob_today)) {
    p_cum[j] <- p_cum[j-1] - ((1 - p_cum[j-1]) * prob_today[j])
  }
  p_cum
}

这仍然有问题。对于j=1，您尝试访问p_cum[0]，它是一个零长度向量，并且您的计算假定一个长度向量。这就是您收到错误消息的原因

Error in p_cum[j] <- p_cum[j - 1] - ((1 - p_cum[j - 1]) * prob_today[j]) : 
  replacement has length zero

初始化p_cum[1]，然后遍历其余部分。

prob_cum <- function(prob_today) {
  p_cum <- rep(0, length(prob_today))
  p_cum[1] <- prob_today[1]
  for (j in 2:length(prob_today)) {
    p_cum[j] <- p_cum[j-1] - ((1 - p_cum[j-1]) * prob_today[j])
  }
  p_cum
}

这种循环结构有潜在的危险。只要prob_today 的长度至少为 2，它就可以工作，但如果长度为 1，它的行为就会出乎意料。更好的是

prob_cum <- function(prob_today) {
  p_cum <- rep(0, length(prob_today))
  p_cum[1] <- prob_today[1]
  for (j in seq_along(prob_today)[-1]) {
    p_cum[j] <- p_cum[j-1] - ((1 - p_cum[j-1]) * prob_today[j])
  }
  p_cum
}

现在我们遇到了一个真正的问题：您的算法是错误的。每天获得至少一场胜利的概率j 是每天获得至少一场胜利的概率j-1 加上当天获得胜利的概率j 假设当时还没有获胜.你有一个减号。

prob_cum <- function(prob_today) {
  p_cum <- rep(0, length(prob_today))
  p_cum[1] <- prob_today[1]
  for (j in seq_along(prob_today)[-1]) {
    p_cum[j] <- p_cum[j-1] + ((1 - p_cum[j-1]) * prob_today[j])
  }
  p_cum
}

现在你有了一个可以工作的函数：

> prob_cum(prob_daily)
[1] 0.500 0.750 0.875
> prob_cum(c(0.5, 0.01, 0.99))
[1] 0.50000 0.50500 0.99505

完全矢量化的解决方案来自于以不同的方式表达概率。至少获得一场胜利的概率是 1 减去到那天为止所有失败的概率。这些是独立的概率，所以只是每天亏损的产物。

prob_cum <- function(prob_today) {
  1 - cumprod(1-prob_today)
}

结果相同

> prob_cum(prob_daily)
[1] 0.500 0.750 0.875
> prob_cum(c(0.5, 0.01, 0.99))
[1] 0.50000 0.50500 0.99505

适用于单个值和空向量，无需任何额外调整

> prob_cum(c(0.75))
[1] 0.75
> prob_cum(c())
numeric(0)