三个或更多集合的可能性

Question

考虑以下几组概率(三个事件不相互排斥):

• 0.05625成功,0.94375失败
• 0.05625成功,0.94375失败
• 0.05625成功,0.94375失败

如何计算至少一个事件发生的概率(即工会)？

如果可能的话,我更喜欢一个通用的,独立的解决方案,它也可以处理4个或更多事件.在这种情况下,我正在寻找的答案是:

0.05625 + 0.05625 + 0.05625 -
0.05625*0.05625 - 0.05625*0.05625 - 0.05625*0.05625 +
0.05625*0.05625*0.05625
##[1] 0.1594358

我的问题最终比标题更广泛,因为我正在寻找可以计算union,intersection(0.05625*0.05625*0.05625 = 0.0001779785),没有事件发生(1 - 0.1594358 = 0.8405642)或恰好发生一个事件(0.150300)的概率的函数.换句话说,这个在线连接三个事件计算器的R解决方案.我已经查看了这个prob包,但是对于这种简单的用例来说,它似乎有一个太复杂的界面.

Answer 1

平等概率

您可以使用二项式密度函数获得正好0,1,2或3 dbinom的概率,二项式密度函数返回在给定独立尝试总数的情况下获得指定成功次数(第一个参数)的概率(第二个参数) )以及每次尝试成功的概率(第三个参数):

dbinom(0:3, 3, 0.05625)
# [1] 0.8405642090 0.1502995605 0.0089582520 0.0001779785

所以,如果你想要至少发生一次的概率,那就是:

sum(dbinom(1:3, 3, 0.05625))
# [1] 0.1594358

要么

1 - dbinom(0, 3, 0.05625)
# [1] 0.1594358

该dbinom功能也可以解决您的其他问题.例如,所有发生的可能性是:

dbinom(3, 3, 0.05625)
# [1] 0.0001779785

正好一个的概率是:

dbinom(1, 3, 0.05625)
# [1] 0.1502996

没有概率是:

dbinom(0, 3, 0.05625)
# [1] 0.8405642

不等概率 - 一些简单的案例

如果在向量中存储了不相等的概率p并且每个项目都是独立选择的,则需要做更多的工作,因为该dbinom函数不适用.仍然,一些计算非常简单.

没有选择项目的概率只是1减去概率的乘积(选择至少一个的概率只是一减去这个):

p <- c(0.1, 0.2, 0.3)
prod(1-p)
# [1] 0.504

所有概率都是概率的乘积:

prod(p)
# [1] 0.006

最后,选择恰好一个的概率是其概率乘以所有其他元素未被选中概率的所有元素的总和:

sum(p * (prod(1-p) / (1-p)))
# [1] 0.398

同样,完全n-1被选中的概率(n概率的数量在哪里)是:

sum((1-p) * (prod(p) / p))
# [1] 0.092

不等概率 - 完整案例

如果你想要每一个可能的成功计数的概率,一个选项可能是计算所有2^n事件组合(这是A. Webb在答案中所做的).相反,以下是O(n ^ 2)方案:

cp.quadratic <- function(p) {
  P <- matrix(0, nrow=length(p), ncol=length(p))
  P[1,] <- rev(cumsum(rev(p * prod(1-p) / (1-p))))
  for (i in seq(2, length(p))) {
    P[i,] <- c(rev(cumsum(rev(head(p, -1) / (1-head(p, -1)) * tail(P[i-1,], -1)))), 0)
  }
  c(prod(1-p), P[,1])
}
cp.quadratic(c(0.1, 0.2, 0.3))
# [1] 0.504 0.398 0.092 0.006

基本上,我们将P_ij定义为我们具有完全i成功的概率,所有这些概率都处于适当位置j或更大.基本情况i=0和i=1计算相对简单,然后我们有以下重复:

P_ij = P_i(j+1) + p_j / (1-p_j) * P_(i-1)(j+1)

在函数中cp.quadratic,我们循环增加i,填充P矩阵(即nx n).因此总操作计数为O(n ^ 2).

例如,这使您能够在一秒钟内计算大量选项的分布:

system.time(cp.quadratic(sample(c(.1, .2, .3), 100, replace=T)))
#    user  system elapsed 
#   0.005   0.000   0.006 
system.time(cp.quadratic(sample(c(.1, .2, .3), 1000, replace=T)))
#    user  system elapsed 
#   0.165   0.043   0.208 
system.time(cp.quadratic(sample(c(.1, .2, .3), 10000, replace=T)))
#    user  system elapsed 
#  12.721   3.161  16.567 

我们可以在一分之一秒内计算1,000个元素的分布,在一分钟内计算10,000个元素的分布; 计算2 ^ 1000或2 ^ 10000个可能的结果将花费相当长的时间(子集的数量分别是301位和3010位数).