Sur un vecteur aléatoire dont les trois composantes suivent des lois binomiales

A, B, C étant trois épreuves de Bernoulli de probabilités de réussite respectives p, q, r (toutes dans ]0;1[), et n étant un entier ≥1, on considére l'épreuve aléatoire E suivante : Le résultat associé à cette réalisation de l'épreuve E est le triplet (ou vecteur) (i,j,k).
On notera Ω(E) l'ensemble des résultats possibles de cette épreuve E.

Sur Ω(E) on définit alors quatre variables aléatoires I, J, K, V qui sont les quatre applications de Ω(E) dans N ou dans N3 suivantes :

On a alors les quatre résultats suivants :

Applications (et origine) de cette étude

Application 1.

Récemment (fin 2011), une entreprise, après modélisation de son problème, a considéré l'épreuve aléatoire E ci-dessus, son but étant de déterminer, à p, q, r fixés, le plus petit entier n tel que P(K≥1)≥90%.
Dans un premier temps elle a utilisé la méthode approchée (trop) décrite à
l'annexe 2 ci-dessous.

En fait, on déduit tout de suite du 3) ci-dessus, que

P(K≥1)=1-(1-pqr)n
et donc
à p, q, r fixés, le plus petit entier n tel que P(K≥1)≥s (avec s dans [0;1[) est le plus petit entier n tel que n≥ln(1-s)/ln(1-pqr).
Deux cas sont à envisager :
Exemples

Remarque : vérification de la formule P(K≥1)=1-(1-pqr)n.

Application 2.

La réalisation pratique des épreuves aléatoires A,B,C ayant des coûts respectifs CA, CB, CC, l'entreprise a cherché à savoir combien allait lui coûter "la" réalisation de E.
Ce coût n'est évidemment pas le même à chaque réalisation de E :
si (i,j,k) est le résultat obtenu lors d'une réalisation de E, c'est qu'il y a eu, outre n réalisations de A, i réalisations de B et j réalisations de C et le coût total de cette réalisation de E est nCA+iCB+ jCC
, qui ne dépend évidemment pas de k.
Ce coût varie entre deux valeurs extrêmes, selon le résultat de E obtenu :


En fait ce coût est la variable aléatoire C définie sur Ω(E) par C=nCA+ICB+ JCC.
cad C(i,j,k)=nCA+iCB+ jCC.

Sachant que l'entreprise cherche à réaliser une seule fois E, avec n assez grand pour lui assurer une forte probabilité que K prenne une valeur ≥1, quel coût doit-elle prévoir?
Le plus réaliste est de prévoir que ce coût sera l'espérance mathématique ou coût moyen de la variable aléatoire C, quantité notée E(C).
C'est la somme de tous les coûts possibles, chaque coût étant pondéré par sa probabilité d'obtention :
E(C)=∑i,j,k/ 0≤k≤j≤i≤nP(i,j,k)C(i,j,k).
Mais l'opérateur espérance E étant linéaire, notamment l'espérance d'une somme est la somme des espérances, on obtient
E(C)=E(nCA+ICB+ JCC)=nCA+E(I)CB+ E(J)CC, soit, en utilisant le 3) ci-dessus

E(C)=n(CA+pCB+pqCC)
Ce coût moyen (ou espérance du coût) est proportionnel à n.

Par exemple, dans le cas p=0.8, q=0.4, r=0.8, si on passe de n=8 à n=9, Compte tenu de ces éléments, c'est alors à l'entreprise de choisir entre n= 8 et n=9.
Preuves

preuve de 1

Si (i,j,k) est un résultat possible de E, nécessairement n≥i≥j≥k≥0 ; mais réciproquement, tout triplet (i,j,k) tel que n≥i≥j≥k≥0 est évidemment un résultat possible de E.
Il s'agit donc de déterminer le nombre N de triplets (ou vecteurs à 3 composantes) (i,j,k) avec n≥i≥j≥k≥0.

1ière méthode ("directe")
Pour cela, on partitionne ces triplets en trois catégories : soit i=j=k, soit i=j≠k ou i≠j=k, soit i>j>k.

2ième méthode : on va se ramener à un résultat classique de dénombrement.
Soit f l'application qui au triplet (i,j,k) avec n≥i≥j≥k≥0 associe le quadruplet (y1, y2, y3, y4) avec


f est évidemment injective.
On remarque que les yi sont des entiers tous ≥0 et leur somme est n ; et étant donné un tel quadruplet, cad (y1, y2, y3, y4) avec yi≥0 et ∑yi=n, il existe un triplet (i,j,k) vérifiant n≥i≥j≥k≥0 et dont l'image par f est justemment ce quadruplet.
En effet, pour cela on prend : Et par sommation de ces trois égalité k=n-y1-y2-y3=y4≥0.
f est donc une bijection de l'ensemble des triplets (i,j,k) avec n≥i≥j≥k≥0 vers l'ensemble des quadruplets d'entiers positifs ou nuls (y1, y2, y3, y4) tels que ∑iyi=n.
Le nombre N de triplets (i,j,k) tels que n≥i≥j≥k≥0 est donc le nombre de quadruplets d'entiers positifs ou nuls (y1, y2, y3, y4) tels que ∑iyi=n : cf un résultat classique il y en a Cn+4-14-1=Cn+33.

preuve de 2

La formule des probabilités conditionnnelles P(A/B)P(B)=P(A ∩ B) donne
P(i,j,k)=P(I=i ∩ J=j ∩ K=k)=P(I=i)P(J=j/I=i)P(K=k/(I=i ∩ J=j)
Or d'où P(i,j,k)=Cni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×Cjk×rk×(1-r)j-k

Remarque : on peut vérifier sans difficulté que la somme S des N probabilités P(i,j,k) est bien 1, cela à partir du fait que i / 0≤i≤mCmi×si×(1-s)n-i=(s+(1-s))m=1.


Cf les six premiers facteurs de P(i,j,k) ne dépendent pas de k cette somme est S=∑i,j/ 0≤j≤i≤nCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×(∑0≤k≤jCjk×rk×(1-r)j-k)
S=∑i,j/ 0≤j≤i≤nCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j
S=∑i/ 0≤i≤nCni×pi×(1-p)n-i×(∑j / 0≤j≤iCij×qj×(1-q)i-j)
S=∑i / 0≤i≤nCni×pi×(1-p)n-i=1

preuve de 3
Vérifions que I suit une loi binomiale de paramètres n et p (puisque ce résultat est évident cf la définition de I).

Il s'agit de vérifier que pour tout i dans {0;1;...;n} on a P(I=i)=Cnipi(1-p)n-i.
P(I=i)=∑j,k / 0≤k≤j≤iCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×Cjk×rk×(1-r)j-k
P(I=i)=Cni×pi×(1-p)n-i×∑j,k / 0≤k≤j≤iCij×qj×(1-q)i-j×Cjk×rk×(1-r)j-k
P(I=i)=Cni×pi×(1-p)n-i×∑j / 0≤j≤iCij×qj×(1-q)i-j×(∑k / 0≤k≤jCjk×rk×(1-r)j-k)
P(I=i)=Cni×pi×(1-p)n-i×∑j / 0≤j≤iCij×qj×(1-q)i-j
P(I=i)=Cni×pi×(1-p)n-i

Montrons que J suit une loi binomiale de paramètres n et pq : là c'est moins évident. Il s'agit de montrer que pour tout j dans {0;1;...;n} on a P(J=j)=Cnj(pq)j(1-pq)n-j.
P(J=j)=∑i,k / 0≤k≤j≤i≤nCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×Cjk×rk×(1-r)j-k
P(J=j)=∑i / j≤i≤nCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×∑k/ 0≤k≤jCjk×rk×(1-r)j-k
P(J=j)=∑i / j≤i≤nCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j
P(J=j)=pjqji / 0≤i-j≤n-jCniCij(p(1-q))i-j(1-p)n-j-(i-j)
Mais CniCij=CnjCn-ji-j
D'où en sortant le Cnj du ∑ ci-dessus on obtient
P(J=j)=pjqjCnji / 0≤i-j≤n-jCn-ji-j(p(1-q))i-j(1-p)n-j-(i-j)
P(J=j)=pjqjCnj(p(1-q)+(1-p))n-j, ce qui donne
P(J=j)=Cnj(pq)j(1-pq)n-j

Montrons que K suit une loi binomiale de paramètres n et pqr : là aussi c'est moins évident. Il s'agit de montrer que pour tout k dans {0;1;...;n} on a P(K=k)=Cnk(pqr)k(1-pqr)n-k.
P(K=k)=∑i,j / k≤j≤i≤nCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×Cjk×rk×(1-r)j-k
P(K=k)=rki / k≤i≤nCni×pi×(1-p)n-i×∑j / k≤j≤iCij×qj×(1-q)i-j×Cjk×(1-r)j-k
Mais CijCjk=CikCi-kj-k=CikCi-ki-j, d'où
P(K=k)=rkqki / k≤i≤nCni×pi×(1-p)n-i×Cikj / 0≤j-k≤i-kCi-kj-k×(q(1-r))j-k×(1-q)i-k-(j-k)
P(K=k)=rkqki / k≤i≤nCni×pi×(1-p)n-i×Cik(q(1-r)+(1-q))i-k
Et CniCik=CnkCn-ki-k donne
P(K=k)=rkqkCnki / k≤i≤nCn-ki-k×pi×(1-p)n-i×(1-qr)i-k
P(K=k)=pkrkqkCnki / 0≤i-k≤n-kCn-ki-k×(p(1-qr))i-k×(1-p)n-k-(i-k)
P(K=k)=pkrkqkCnk(p(1-qr)+(1-p))n-k, ce qui donne
P(K=k)=Cnk(pqr)k(1-pqr)n-k

preuve de 4

P(I=n ∩ J=n ∩ K=n)=P(n,n,n)=(pqr)n≠P(I=n)P(J=n)P(K=n)=pn(pq)n(pqr)n : donc les variables aléatoires I, J, K ne sont pas indépendantes mutuellement. On verra en outre, ci-dessous, qu'elles ne sont pas indépendantes 2 à 2.

Cf les espérances de I, J, K sont connues, il y a juste à calculer celles de IJ, IK, JK

On sera amener à utiliser les résultats suivants : si X suit une loi binomiale de paramètres m et s, alors


Covariance de I et J E(IJ)=∑i,j,k / 0≤k≤j≤i≤nijP(i,j,k)
E(IJ)=∑i,j / 0≤j≤i≤nijCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×∑k / 0≤k≤jCjk×rk×(1-r)j-k
E(IJ)=∑i,j / 0≤j≤i≤nijCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j
E(IJ)=∑i / 0≤i≤niCni×pi×(1-p)n-i×∑j / 0≤j≤ijCij×qj×(1-q)i-j
E(IJ)=∑i / 0≤i≤niCni×pi×(1-p)n-i×(iq), d'après R1 et donc
E(IJ)=qE(I2), ce qui donne cf R2
E(IJ)=q(np(1-p)+n2p2) et
cov(I,J)=E(IJ)-E(I)E(J)=q(np-np2+n2q2)-np×npq=npq(1-p)
Cette covariance étant non nulle, I et J ne sont pas indépendantes
Covariance de I et K E(IK)=∑i,j,k / 0≤k≤j≤i≤nikP(i,j,k)
E(IK)=∑i,j / 0≤j≤i≤niCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×∑k / 0≤k≤jkCjk×rk×(1-r)j-k
E(IK)=∑i,j / 0≤j≤i≤niCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×(jr) d'après R1 et donc
E(IK)=r∑i / 0≤i≤niCni×pi×(1-p)n-i×∑j / 0≤j≤ijCij×qj×(1-q)i-j
E(IK)=r∑i / 0≤i≤niCni×pi×(1-p)n-i×(iq) d'après R1 et
E(IK)=rq∑i / 0≤i≤ni2Cni×pi×(1-p)n-i
E(IK)=rqE(I2)=rq(np(1-p)+n2p2)
cov(I,K)=E(IK)-E(I)E(K)=rq(np-np2+n2q2)-np×npqr=npqr(1-p)
Cette covariance étant non nulle, I et K ne sont pas indépendantes
Covariance de J et K E(JK)=∑i,j,k / 0≤k≤j≤i≤njkP(i,j,k)
E(JK)=∑i,j / 0≤j≤i≤njCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×∑k / 0≤k≤jkCjk×rk×(1-r)j-k
E(JK)=∑i,j / 0≤j≤i≤njCni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×(jr) d'après R1 et
E(JK)=r∑i / 0≤i≤nCni×pi×(1-p)n-i×∑j / 0≤j≤ij2Cij×qj×(1-q)i-j
E(JK)=r∑i / 0≤i≤nCni×pi×(1-p)n-i×(iq(1-q)+i2q2) d'après R2 et
E(JK)=r[q(1-q)∑i / 0≤i≤niCni×pi×(1-p)n-i+q2i / 0≤i≤ni2Cni×pi×(1-p)n-i]
E(JK)=r[q(1-q)np+q2(np(1-p)+n2p2]
E(JK)=r(npq(1-pq)+(npq)2) ; c'est rE(J2) !
cov(J,K)=E(JK)-E(J)E(K)=r(npq(1-pq)+(npq)2)-npq×npqr=npqr(1-pq)
Cette covariance étant non nulle, J et K ne sont pas indépendantes
Matrice de covariance de V, V étant considérée comme matrice ligne

La définition donnée en remarque 2 de l'énoncé permet de dire tout de suite que cette matrice de covariance est la matrice Δ symétrique suivante (écrite sous forme de tableau) :

np(1-p)npq(1-p)npqr(1-p)
npq(1-p)npq(1-pq)npqr(1-pq)
npqr(1-p)npqr(1-pq)npqr(1-pqr)

Elle est factorisable par np.

Montrons que Δ est définie positive.
Soit x une matrice colonne à trois éléments, tous réels : x est une "constante".
txΔx=E(tx[t(V-E(V))(V-E(V))]x)
txΔx=E([txt(V-E(V))][(V-E(V))x])
txΔx=E(z2), z étant la variable aléatoire réelle txt(V-E(V))=(V-E(V))x.
Donc pour tout x, txΔx≥0 : Δ est positive.
Il faut montrer maintenant que si x≠0, txΔx>0.
Suposons que txΔx=0, cad que E(z2)=0 :
cela équivaut à z est la variable aléatoire ne pouvant prendre qu'une valeur avec une probabilité non nulle : c'est la valeur 0, prise avec la probabilité 1.
En effet si z prend une valeur non nulle avec une probabilité non nulle, alors E(z)>0.
En notant tx=(x1 x2 x3), z=x1(I-E(I))+x2(J-E(J)) +x3(K-E(K)),
soit z=x1I+x2J +x3K-r avec r=x1E(I)+x2E(J) +x3E(K).
D'où


si V prend la valeur (1,1,1), z prend la valeur x1+x2+x3-r
si V prend la valeur (1,1,0), z prend la valeur x1+x2-r
si V prend la valeur (1,0,0), z prend la valeur x1-r
si V prend la valeur (0,0,0), z prend la valeur -r
Or, cf la remarque du 2), la probabilité que z prenne chacune de ces quatre valeurs est non nulle, donc z ne pouvant prendre que la valeur 0 avec une probabilité non nulle, c'est que r=0, puis x1=x2=x3=0, soit x=0.
Donc on a montré txΔx=0 => x=0, cad x≠0 => txΔx>0 : Δ est bien définie positive.

Pour calculer son déterminant (qui doit être >0 puisque cette matrice de covariance est définie-positive), on commence évidemment par la mise en facteur de np, puis on factorise 1-p dans la 1ière ligne, q dans la seconde et qr dans la 3ième.
Ainsi det(Δ)=(np)3q2r(1-p)det(Δ') avec Δ' la matrice ci-dessous

1qqr
1-p1-pqr(1-pq)
1-p1-pq1-pqr

En ajoutant aux 2ième et 3ième lignes la 1ière multipliée par p, on obtient que det(Δ') est le déterminant de la matrice ci-dessous

1qqr
11r
111
Sarrus donne alors tout de suite det(Δ')=1+2qr-qr-r-q=1-r-q+qr=(1-q)(1-r).
Finalement le déteminant de Δ est det(Δ)=n3p3q2r(1-p)(1-q)(1-r)>0.

Quant aux valeurs propres de la matrice de covariance Δ, ce sont les valeurs propres de Δ''=(1/np)Δ, multipliées par np.
Comme Δ'' ne dépend pas de n, les valeurs propres de Δ'' ne dépendent pas aussi de n, donc les valeurs propres de Δ sont proportionnelles à n.

Remarque : dans le cas p=q=r=0.5, la matrice Δ'' est

1/21/41/8
1/43/83/16
1/83/167/32

Son polynôme caractéristique est -X3+(35/32)X2-(17/64)X+1/64.
Δ'', donc Δ, n'admettent aucune valeur propre rationnelle car si a/b est racine de ce polynôme caractéristique, avec a et b dans Z et premiers entre eux, alors -64a3+70a2b-17ab2+b3=0 ; donc a doit diviser b, soit |a|=1, et b doit diviser 64.
Les seules possibilités pour une racine rationnelle sont alors ±2-i pour i=0, 1, 2, 3, 4, 5, 6 : aucune n'est effectivement racine.
En fait les valeurs propres de Δ'' sont approximativement 0.77823, 0.22712, 0.0884 et donc les valeurs propres de la matrice de covariance Δ sont approximativement (n/2)×0.77823, (n/2)×0.22712, (n/2)×0.0884.
Leur somme est, approximativement (n/2)×1.09375, ce qui était attendu car la trace de la matrice de covariance, qui est la somme de ses termes diagonaux, cad la somme des variances de I, J, K est, np(1-p+q(1-pq)+qr(1-pqr))=(n/2)×(35/32)=(n/2)×1.09375.
Et leur produit est approximativement 0.0156248(n/2)3, ce qui là aussi était attendu, puisque le déterminant de Δ est n3/29=0.015625(n/2)3.

preuve de 5

φV(T)=i,j / 0≤j≤i≤nexp(it1i)exp(it2j)Cni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×∑k / 0≤k≤jCjk×(rexp(it3))k×(1-r)j-k
Et par utilisations successives de la formule du binôme,
φV(T)=∑i,j / 0≤j≤i≤nexp(it1i)exp(it2j)Cni×pi×(1-p)n-i×Cij×qj×(1-q)i-j×(1-r+rexp(it3))j
φV(T)=∑i / 0≤i≤nexp(it1i)Cni×pi×(1-p)n-i×∑j / 0≤j≤iexp(it2j)Cij×qj×(1-q)i-j×(1-r+rexp(it3))j
φV(T)=∑i / 0≤i≤nexp(it1i)Cni×pi×(1-p)n-i×∑j / 0≤j≤iCij×((1-r)exp(it2)+rexp(i(t2+t3)))jqj×(1-q)i-j
φV(T)=∑i / 0≤i≤nexp(it1i)Cni×pi×(1-p)n-i×(1-q+q(1-r)exp(it2)+qrexp(i(t2+t3)))i
φV(T)=∑i / 0≤i≤nCni×pi×(1-p)n-i×((1-q)exp(it1)+q(1-r)exp(i(t1+t2))+qrexp(i(t1+t2+t3)))i
φV(T)=∑i / 0≤i≤nCni×(p(1-q)exp(it1)+pq(1-r)exp(i(t1+t2))+pqrexp(i(t1+t2+t3)))i(1-p)n-i
φV(T)=[1-p+p(1-q)exp(it1)+pq(1-r)exp(i(t1+t2))+pqrexp(i(t1+t2+t3))]n

On vérifie


Annexe 1

Calcul numérique en ligne de P(K≥1)

N'ayant pas vu tout de suite, comme "l'entreprise", que P(K≥1) est en fait égal à 1-(1-pqr)n, je m'étais "amusé" à programmer un calcul numérique en ligne de P(K≥1) à partir du fait que cette probabilité est évidemment la somme des P(i,j,k), sommation prise sur les (i,j,k) tels que n≥i≥j≥k≥1.
Comme pratiquement il n'y a que moi qui me lit, je garde ce module en souvenir ...

Rentrer les quatre données dans les champs ci-dessous :

<== mettre dans ce champ la valeur de n≥1

<== mettre dans ce champ la valeur de p (un point pour la virgule : par exemple 0.8)

<== mettre dans ce champ la valeur de q

<== mettre dans ce champ la valeur de r

et lancez le calcul en cliquant sur le bouton ...lancement du calcul de P(K≥1)!

Le calcul de la somme des P(i,j,k) concernés est fait avec un petit programme en Java-Script, lequel doit donc être activé dans votre navigateur.

Valeur de P(K≥1)

Et bien entendu rien n'empêche de rentrer un autre jeu de données (on peut ne modifier qu'une des données).

Application : ce petit module de calcul permet d'obtenir très rapidement, en quelques clics, le plus petit entier n ( à p,q,r donnés) tels que P(K≥1)≥0.9 (par exemple) :

Note :
ce module de calcul a été testé dans les cas suivants (ne pas hésiter à m'écrire en cas de problème : adresse sur la page d'accueil) :


Annexe 2retour

Sur une approximation du plus petit entier n ( à p, q, r fixés ) tel que lorsqu'on réalise l'épreuve ci-dessus, on ait P(K≥1)≥0.9 (par exemple)

Ce plus petit entier n ( à p, q, r fixés ) tel que lorsqu'on réalise l'épreuve E ci-dessus, on ait P(K≥1)≥0.9 sera noté nmin.

Voici la méthode pour approximer nmin :

C'est cet entier n' que l'entreprise envisageait de prendre comme approximation de nmin.

Avant de juger de la qualité de cette approximation, précisons la méthode pratique pour obtenir cet entier n'.

Soit Xm,s une variable aléatoire suivant une loi binomiale de paramètres m et s, et Fm,s sa fonction de répartition : Fm,s(k)=P(Xm,s≤k) et Fm,s(0)=(1-s)m, Fm,s(m)=1.
Puisque, pour k entier, P(Xm,p≥k)=1-Fm,p(k-1),

Comme pratiquement dans tout ouvrage sur les variables aléatoires réelles on trouve des tables de Fm,s, ces dernières peuvent être utilisées pour déterminer i' et n' : mais bien souvent ces tables ne sont pas par pas de 1 en m, et les valeurs de s proposées ne sont pas forcément celles auxquelles on s'intéresse (p, q, r), et donc il faut faire des interpollations : il vaut mieux aller sur calcul en ligne de la loi binomiale où il y a un module de calcul qui après avoir rentré n'importe quel m et et n'importe quel s donne "toute" la fonction de répartition Fm,s, cad il donne les valeurs de Fm,s(k), pour k=0 à m : j'ai utilisé ce lien pour faire les exemples ci-dessous.

Sur deux exemples comparons cet entier n' avec nmin.

Exemple 1 : p=q=r=0.5

alors que dans ce cas nmin=18, d'où une erreur relative |n'-nmin|/nmin≈66%.

Exemple 2 : p=0.8, q=0.4, r=0.8

alors que dans ce cas nmin=8, d'où une erreur relative |n'-nmin|/nmin≈62%.

Conclusion : cette approche pour approximer nmin est mauvaise car elle peut conduire à une erreur relative de plus de 60%.

La raison est que cette méthode approchée s'impose des conditions qui ne sont pas nécessaires.
Par exemple lorsque p=0.8, q=0.4, r=0.8 cette méthode approchée cherche le plus petit entier n'≥9 tel qu'en réalisant n' fois l'épreuve A, la probabilité d'avoir au moins 9 réussites soit ≥0.9 : or l'étude de l'épreuve aléatoire E, montre que pour n=8 on a P(K≥1)≥0.9, sans pour autant avoir P(I≥9)≥0.9, puisque P(I≥9)=0, la valeur maximum de I étant alors 8!