按结果分组的事件概率

Probabilities of events grouped by their outcome

假设有 $n$ 个独立事件。每个都有概率 $p_n$ 和相关损失 $l_n$。我的目标是列出所有可能的损失金额及其相关概率。

最终我想将其扩展到具有可变概率和损失金额的 10-20 个事件集。这一切都将在 R 中完成。

各种结果由幂集给出,例如对于三个事件:(null)、(A)、(B)、(C)、(A 和 B)、(A 和 C)、(B 和 C)、(A 和 B 和 C)。这些结果中的每一个的概率可以通过取每个子集中的概率的乘积来找到,总损失可以通过取每个子集中的损失总和来找到。

我的问题是如何按损失量进行汇总,即找到幂集中所有唯一的损失量并产生它们的概率。

我觉得 inclusion/exclusion principle 已经完成了一半,但我不太清楚如何将它应用于我的特定问题,尤其是当事件数量超过 3 时,或者对于中等大小的集合,例如如何对以上所有 2 个元素集进行分组。

对于这么小的问题——最多有 2^20(大约一百万)种可能性——蛮力就可以了。

为了说明,让我们生成一些中等大小的数据:

n <- 15
set.seed(17)
p <- runif(n)
loss <- ceiling(rgamma(n, 3, 1/2))
signif(rbind(Probability=p, Loss=loss), 2)

以下是此示例的输入值:

Probability  0.16 0.97  0.47 0.78  0.41 0.54  0.21 0.19 0.78  0.19  0.43 0.0023  0.83  0.83  0.96
Loss        12.00 4.00 10.00 8.00 10.00 6.00 12.00 5.00 4.00  8.00  8.00 8.0000  4.00  4.00  4.00

expand.grid生成幂集的二进制指标,然后使用数组运算相对快速地计算损失和所有可能结果的概率:

powerset <- t(expand.grid(lapply(p, function(x) 0:1)))
probability <- apply(powerset * (2*p - 1) + (1-p), 2, prod)
losses <- colSums(powerset * loss)

(在这个老化的 Xeon 工作站上,当 n 为 20 时,这最多需要 5 秒。)

使用 tapply 按损失汇总:

x <- tapply(probability, losses, sum)

(当 n 为 20 时,这又需要 1 到 2 秒。)

我们可以通过 (a) 验证概率之和是否为 1 以及 (b) 检查预期损失是否为各个事件的预期损失之和来检查一致性:

if(sum(probability) - 1 != 0) warning("Unnormalized probability.")
if(sum(probability * losses) - sum(p*loss) != 0) warning("Inconsistent result.")

让我们绘制最终的损失分布图。

library(ggplot2)
ggplot(data.frame(Loss=as.numeric(names(x)), Probability=x), 
       aes(Loss, Probability)) + 
  geom_col(color="White")