Probabilités d'événements regroupés selon leur résultat

Question

Disons qu'il y a $n$ événements indépendants. Chacun a une probabilité $p_n$ et une perte associée $l_n$. Mon objectif est de produire une liste de tous les montants de perte possibles et de leurs probabilités associées.

À terme, j'aimerais étendre ce système à des ensembles de 10 à 20 événements avec des probabilités et des montants de pertes variables. Tout cela sera fait en R.

Les différentes issues sont données par l'ensemble de puissance, par exemple pour trois événements : (nul), (A), (B), (C), (A et B), (A et C), (B et C), (A et B et C). La probabilité de chacune de ces issues peut être trouvée en prenant le produit des probabilités dans chaque sous-ensemble, et la perte totale en prenant la somme des pertes dans chaque sous-ensemble.

Mon problème est de savoir comment agréger les montants de perte, c'est-à-dire trouver tous les montants de perte uniques dans l'ensemble de puissance et produire leurs probabilités.

J'ai l'impression d'avoir fait la moitié du chemin avec les principe d'inclusion/exclusion Mais je n'arrive pas à comprendre comment l'appliquer à mon problème particulier, surtout lorsque le nombre d'événements dépasse 3, ou dans le cas des ensembles de taille intermédiaire, par exemple comment regrouper tous les ensembles à 2 éléments ci-dessus.

Answer 1

Pour un problème aussi petit - il y a au maximum 2^20 (environ un million) de possibilités - la force brute fonctionne bien.

Pour illustrer, générons quelques données de taille moyenne :

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

Voici les valeurs d'entrée pour cet exemple :

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

Générer un indicateur binaire de l'ensemble de puissance avec expand.grid puis utiliser les opérations de tableau pour un calcul relativement rapide des pertes et des probabilités de tous les résultats possibles :

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

(Sur cette station de travail Xeon vieillissante, cela prend jusqu'à 5 secondes lorsque n est de 20).

Résumez par perte en utilisant tapply :

x <- tapply(probability, losses, sum)

(Cela prend encore 1 à 2 secondes lorsque n est de 20).

Nous pouvons vérifier la cohérence en (a) vérifiant que la somme des probabilités est égale à l'unité et (b) en vérifiant que la perte attendue est la somme des pertes attendues des événements individuels :

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

Traçons la distribution des pertes qui en résulte.

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