0

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

假设我每天要玩一次游戏，并且我知道每天获胜的概率。我想要一个函数，它采用概率向量并返回至少一天成功的累积概率。因此，如果我连续玩 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)
4

2 回答 2

5 
>  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
于 2013-03-27T20:29:03.577 回答
4

一次解决每个问题：

你有一个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)
于 2013-03-27T20:54:53.740 回答