Convolution de variable aléatoire.

La convolution est un concept qui a une apparence plutôt simple de loin. L’idée est très simple : si j’ai deux variables aléatoires et que je les additionne, qu’est-ce que j’obtiens à la fin ?

Cela semble plutôt simple comme question à première vue, mais lorsque l’on se penche sur l’application et les calculs, on se rend vite compte que la question est loin d’être triviale. Personnellement, je me suis cassé la tête sur ce problème pendant beaucoup plus d’heures que je ne voudrais l’admettre…

Donc l’idée va être de simplifier le concept au maximum tout en l’explorant de façon très complète.

I - Introduction

Insane Monopoly >>> Monopoly :)

Étant donné que le sujet est déjà bien abstrait comme ça, je vais utiliser le Monopoly comme exemple afin de mieux visualiser.

Lorsqu’on joue au Monopoly, n’importe quelle personne correctement constituée dit : « Je dois avancer de 8 cases » et non « Je dois avancer de 3 cases puis de 5 cases ». Pourtant, lorsque l’on modélise mathématiquement le tirage des dés, on utilise deux variables aléatoires distribuées de façon uniforme et discrète, indépendantes l’une de l’autre.

$E = \{1,2,3,4,5,6\}$

$X \sim \mathcal{U}(E) \quad \text{et} \quad Y \sim \mathcal{U}(E)$

$\mathbb{P}(X = k) = \frac{1}{6} \quad \text{et} \quad \mathbb{P}(Y = k) = \frac{1}{6} \quad \text{pour } k \in E$

Bien que l’on joue au Monopoly avec deux dés, ce qui compte à la fin, c’est la somme des deux dés. Et la convolution, c’est exactement cela : additionner deux variables aléatoires.

Alors cela peut paraître évident, mais la probabilité de faire $7$ au Monopoly, par exemple, n’est pas du tout $\frac{1}{12}$, mais le résultat d’un calcul de probabilité obtenu via une convolution.

Il est tout à fait possible de le faire « à l’arrache », sans trop de rigueur, en ayant l’intuition suivante :

« En comptant toutes les paires menant à $7$ et en additionnant leurs probabilités, je peux trouver la probabilité de faire $7$. »

Prenons l’exemple de la probabilité de faire $3$ au Monopoly. Pour faire $3$, je peux soit tirer $(1,2)$, soit $(2,1)$. Chacune des paires a une probabilité de $\frac{1}{36}$, donc la probabilité de faire $3$ est

$\frac{1}{36} + \frac{1}{36} = \frac{1}{18}.$

Ici, on l’a fait un peu « à l’arrache », mais il existe une formule qui condense tout ce calcul en une seule ligne.

$p_Z(z) = \sum_{x \in R_X} p_X(x)\, p_Y(z - x)$

Elle peut sembler un peu chaotique à première vue, mais elle est tout à fait intuitive. Voici son fonctionnement :

Prenons toutes les issues possibles du premier dé : 1,2,3,4,5 et 6.

Pour chaque itération, on multiplie la probabilité de cette issue (via la pmf de $X$ ici) avec la probabilité de l’issue désirée de la convolution moins celle obtenue.

Pour l’exemple de la convolution à 3, si on a pioché 1 dans le premier dé, le seul moyen d’atteindre 3 est de piocher 2. On fait donc $3 - 1 = 2$.

C’est pour cela qu’on met $z - x$ dans la pmf de $Y$. On fait simplement cela pour toutes les issues possibles du premier dé de la première variable aléatoire.

Notez d’ailleurs un comportement intéressant : si on pioche, disons, 5 au premier dé et qu’on veut obtenir 3 au total, dans la pmf de $Y$ on va y mettre $3 - 5$, soit $-2$. Or $-2$ n’est pas une issue possible de la variable aléatoire $Y$, donc sa pmf va retourner 0, et comme on fait un produit, cette partie de la somme sera nulle. Cela permet de ne pas sortir du champ des possibles.

Mais alors comment fait-on si on veut additionner 4 ou $n$ variables aléatoires ? Car avec cette formule on ne peut le faire qu’entre deux variables aléatoires.

Prenons l’exemple avec $X + Y + Z$ :

$X + Y = B$

$B + Z = C$

Il suffit simplement de boucler autant de fois que nécessaire et de faire une convolution sur le résultat d’une convolution précédente.

J'ai par ailleurs visualiser une convolution des des du monopoly via matplotlib sur un script Python.

Alors c’est plutôt flagrant : on obtient une forme bien connue, une sorte de gaussienne ! Et c’est l’une des propriétés fondamentales de la convolution : au bout de n convolutions (sous certaines conditions), on tend vers une forme très proche de la gaussienne, peu importe la distribution initiale des variables aléatoires.

C’est d’ailleurs ce qui explique en partie son omniprésence dans la nature. Prenons la taille des êtres humains : le corps humain étant un système bien complexe, il dépend d’un grand nombre de variables aléatoires de types de distributions variés. La taille peut alors être vue comme « une sorte de convolution » de toutes ces variables. Et comme la convolution a tendance à produire une forme gaussienne lorsque l’on additionne un grand nombre de variables indépendantes, cela explique pourquoi la loi normale apparaît si souvent.

II - Le cas continu

Tout ce que nous avons vu jusque-là n’était qu’un amuse-bouche. C’est à partir d’ici que, personnellement, je suis resté coincé pendant beaucoup plus d’heures que je ne voudrais l’admettre. En effet, la technique consistant à compter toutes les combinaisons ne fonctionne plus ! Dans le cas continu, il existe une infinité de combinaisons possibles. Et lorsqu’on touche à l’infini, l’intégrale n’est jamais bien loin.

Sauf que, comme nous allons le voir dans quelques instants, cette intégrale peut être particulièrement chiante à calculer. Et la cerise sur le tas de m#rde, c’est que les bornes sont elles aussi délicates à déterminer, car il peut y avoir plusieurs intégrales selon la zone dans laquelle on se trouve !

Fort heureusement, je vais vous présenter un algorithme très mécanique qui vous permettra toujours de retomber sur vos pattes.

Commençons par un exemple. Imaginons un nouveau scénario :

Nous sommes une start-up de maccha-late ( saveurs fruit du dragon, ice-kiwi je ne sais trop quoi) pour bobos parisiens et ados sans personnalités et nous devons modéliser la quantité de graines de café que nos fournisseurs risquent de nous fournir. Nous avons deux variables aléatoires $C$ et $V$ représentant le nombre de tonnes de graines de café que nous allons recevoir sur l’année, provenant respectivement de la Colombie et du Vietnam. Bien sûr, comme ces deux pays sont très éloignés, nous admettrons l’indépendance entre ces deux variables aléatoires.

$C \sim \mathcal{U}([0,1]) \qquad \text{et} \qquad V \sim \mathcal{U}([0,2])$

Nous cherchons alors à déterminer la loi de

$Z = C + V.$

Ce qui nous intéresse vraiment, c’est ce que nous allons avoir en stock. Nous voulons donc déterminer Z.

Pour se faire nous allons employer la formule de convolution mais continu cette fois.

$F_Z(z) = \mathbb{P}(Z \le z).$

Par définition de $Z$ :

$F_Z(z)=\mathbb{P}(C+V \le z) = \int_{-\infty}^{\infty}\int_{-\infty}^{z-v} f_C(c)\,f_V(v)\,dc\,dv = \int_{-\infty}^{\infty} F_C(z-v)\,f_V(v)\,dv.$

On remarquera qu’ici, j’ai préféré donner la CDF et non la PDF, contrairement au cas discret où nous faisions nos calculs avec la PMF et où il suffisait ensuite de sommer la PMF pour obtenir une probabilité via une inégalité. Si j’ai choisi la CDF, c’est tout simplement parce qu’il est très facile d’en déduire la PDF : la PDF est la dérivée de la CDF par rapport à $z$.

$f_Z(z)=\frac{d}{dz}F_Z(z)$

$f_Z(z)=\frac{d}{dz}\left(\int_{-\infty}^{\infty} F_C(z-v)\,f_V(v)\,dv\right)$

$f_Z(z)=\int_{-\infty}^{\infty} \frac{d}{dz}F_C(z-v)\,f_V(v)\,dv$

$f_Z(z)=\int_{-\infty}^{\infty} f_C(z-v)\,f_V(v)\,dv$

On voit tout de suite le parallèle avec le cas discret. Avançons dans le processus et remplaçons les PDF par leur véritable définition. Les densités de variables aléatoires uniformément distribuées sont relativement simples : ce sont de simples constantes.

Comme

$C \sim \mathcal{U}([0,1]) \qquad \text{et} \qquad V \sim \mathcal{U}([0,2]),$

on a

$f_C(c)= \begin{cases} 1 & \text{si } 0 \le c \le 1 \\ 0 & \text{sinon} \end{cases} \qquad f_V(v)= \begin{cases} \frac{1}{2} & \text{si } 0 \le v \le 2 \\ 0 & \text{sinon} \end{cases}$

Sauf que le problème, c’est que nous avons des fonctions définies par morceaux.

Donc nous aurons également une convolution définie par morceaux. Pour pouvoir nous en sortir, il faut déterminer le domaine de $z$ pour lequel l’intégrande n’est pas nul.

Commençons par déterminer le domaine de

$f_C(z-v)=1 \quad \text{si} \quad 0 \le z-v \le 1.$

Cela donne :

$0 \le z-v \le 1$

$- z \le - v \le 1 - z$

$z-1 \le v \le z.$

Nous connaissons désormais le domaine de $f_C(z-v)$, mais aussi celui de $f_V(v)$ :

$0 \le v \le 2.$

Le domaine d’intégration est donc l’intersection :

$v \in [0,2] \cap [z-1,z].$

(On s'en fou d'integrer ailleurs car c'est tout le temps nulle)

Les bornes deviennent alors :

$\text{borne inférieure} = \max(0,z-1), \qquad \text{borne supérieure} = \min(2,z).$

Etant donner qu'il y'a des fonctions min et max on va devoir trouver leurs points de bacules.

On résout :

$z-1 > 0 \quad \Longleftrightarrow \quad z>1,$

$z < 2.$

On obtient alors trois cas :

$\text{1) } 0 \le z \le 1, \qquad \text{2) } 1 \le z \le 2, \qquad \text{3) } 2 \le z \le 3.$

Il ne nous reste plus qu’à poser les intégrales en faisant varier les bornes selon ces points de bascule.

$f_Z(z)=\frac{1}{2} \int_{\max(0,z-1)}^{\min(2,z)} dv.$

$\textbf{Cas 1 : } 0 \le z \le 1$

$f_Z(z)=\frac{1}{2}\int_{0}^{z} dv =\frac{1}{2}z.$

$\textbf{Cas 2 : } 1 \le z \le 2$

$f_Z(z)=\frac{1}{2}\int_{z-1}^{z} dv =\frac{1}{2}(z-(z-1)) =\frac{1}{2}.$

$\textbf{Cas 3 : } 2 \le z \le 3$

$f_Z(z)=\frac{1}{2}\int_{z-1}^{2} dv =\frac{1}{2}(2-(z-1)) =\frac{3-z}{2}.$

On obtient finalement :

$f_Z(z)= \begin{cases} \dfrac{z}{2} & \text{si } 0 \le z \le 1, \\[8pt] \dfrac{1}{2} & \text{si } 1 \le z \le 2, \\[8pt] \dfrac{3-z}{2} & \text{si } 2 \le z \le 3, \\[8pt] 0 & \text{sinon}. \end{cases}$

La méthode fonctionne avec n’importe quel type de variables continues. J’ai choisi des lois uniformes parce que ce sont celles qui illustrent le plus simplement ce problème de bornes. Mais évidemment, on peut appliquer la méthode à n’importe quel type de variable aléatoire.

Pour mieux comprendre « visuellement » comment fonctionne cette intégration, je ne peux que vous conseiller la superbe vidéo de 3Blue1Brown à ce sujet. (Elle m’a beaucoup aidé à mieux comprendre.)

III - Trucs & Astuces

Comme on a pu le voir précédemment, la convolution est une sacrée horreur à calculer. Le nombre d’erreurs d’inattention potentielles peut être assez affligeant. Fort heureusement, on peut parfois esquiver la méthode « pure ».

Cas numéro 1 — Addition de gaussiennes

On additionne deux variables gaussiennes indépendantes :

$X \sim \mathcal{N}(\mu_X,\sigma_X^2), \qquad Y \sim \mathcal{N}(\mu_Y,\sigma_Y^2).$

Alors,

$Z = X+Y \sim \mathcal{N}(\mu_X+\mu_Y,\ \sigma_X^2+\sigma_Y^2).$

Ici, rien de plus simple : on additionne les moyennes et les variances, et le tour est joué.

Cas numéro 2 — Vous avez un ordinateur sur vous

L’une des méthodes les plus simples, et très répandue grâce à sa simplicité, est la méthode de Monte-Carlo.

La technique est plutôt simple : on génère un grand nombre d’échantillons,

$(X_1,Y_1), \dots, (X_N,Y_N),$

où chaque couple est indépendant et suit les lois de $X$ et $Y$.

On définit alors

$Z_i = X_i + Y_i.$

Pour estimer une probabilité du type

$\mathbb{P}(Z \le z),$

on utilise :

$\mathbb{P}(Z \le z) = \mathbb{E}\!\left[\mathbf{1}_{X+Y \le z}\right]$

et son approximation Monte-Carlo :

$\mathbb{P}(Z \le z) \approx \frac{1}{N} \sum_{i=1}^{N} \mathbf{1}_{X_i+Y_i \le z}.$

La méthode fonctionne grâce à la loi des grands nombres, qui s’applique lorsque l’on répète une expérience un très grand nombre de fois.

Lorsque les bornes ou les pdfs sont une véritable purge, ne cherchez pas plus loin : utilisez Monte-Carlo, cela fonctionne très bien. Personnellement, c’est ma méthode favorite ; je m’embête très rarement à calculer la convolution moi-même.

(C’est un comble qu’un post sur la convolution se termine par une incitation à ne pas utiliser la convolution.)