r - 三组或多组并集的概率

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

我的问题最终比标题更广泛，因为我正在寻找可以计算并集、交集( 0.05625*0.05625*0.05625 = 0.0001779785 )、没有事件发生( 1 - 0.1594358 = 0.8405642 )或仅发生一个事件( 0.150300 )概率的函数)。换句话说，此在线 Conjunction of three events calculator 的 R 解决方案.我已经查看了 prob包，但对于这样一个简单的用例来说，它的接口(interface)似乎太复杂了。

Answer 1

等概率

您可以使用二项式密度函数 dbinom 获得其中恰好 0、1、2 或 3 个发生的概率。 ，它返回在给定独立尝试的总数(第二个参数)和每次尝试的成功概率(第三个参数)的情况下准确获得指定成功次数(第一个参数)的概率:

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 减去概率的乘积(至少一个被选中的概率就是 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矩阵(即 n x 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 位的数字)。