三个或更多集合的并集概率

Question

三个或更多集合的并集概率

rprobabilityset-intersectionset-union

7

考虑以下概率集合（三个事件不互斥）：

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）概率的函数。换句话说，我需要一个R解决方案来处理在线三个事件的结合计算器。我已经研究了prob包，但它似乎对于这样一个简单的用例来说界面太复杂了。

- landroni

1

所有事件的概率相同吗？ - Heroka

1

@Heroka 在这个特定的例子中，是的。但它们不必如此。我正在寻找一种通用解决方案，当所有事件的概率可以不同时。 - landroni

2个回答

3

这里有一个函数可以创建所有事件组合，计算它们的概率，并按出现次数进行聚合。

cp <- function(p) 
{
  ev <- do.call(expand.grid,replicate(length(p),0:1,simplify=FALSE))
  pe <- apply(ev,1,function(x) prod(p*(x==1)+(1-p)*(x==0)))
  tapply(pe,rowSums(ev),sum)
}

使用独立发生概率为0.1、0.2和0.3的事件的概率，与josilber的示例相同：

cp(c(0.1,0.2,0.3))

    0     1     2     3 
0.504 0.398 0.092 0.006

例如，恰好发生两个独立事件的概率为0.092。

注：该文本是关于概率论和统计学中的概率分布的简单示例。

- A. Webb

1

因此，使用cp()方法，至少发生一个事件的概率为：sum(cp(c(0.1,0.2,0.3))[-1]) = 0.496。谢谢！ - landroni

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- josliber · Accepted Answer

等概率

您可以使用二项式密度函数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

恰好为1的概率是：

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 时，恰好选中 n-1 的概率为：

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位数）。