Loi de Dirichlet

Un article de Wikipédia, l'encyclopédie libre.
Aller à : navigation, rechercher
Plusieurs images de la densité de la loi de Dirichlet distribution lorsque K=3 pour différents vecteurs de paramètres α. Dans le sens horaire à partir du coin supérieur gauche : α=(6, 2, 2), (3, 7, 5), (6, 2, 6), (2, 3, 4).

En probabilité et statistiques, la loi de Dirichlet, souvent notée Dir(α), est une famille de lois de probabilité continues pour des variables aléatoires multinomiales. Cette loi (ou encore distribution) est paramétrée par le vecteur α de nombres réels positifs et tire son nom de Johann Peter Gustav Lejeune Dirichlet. Elle est vue comme la généralisation multinomiale de la loi bêta.

Densité de probabilité[modifier | modifier le code]

Illustration du changement de la log-densité avec K=3 lorsque le vecteur α varie de α=(0.3, 0.3, 0.3) à (2.0, 2.0, 2.0), en conservant toutes les composantes individuelles \alpha_i égales entre elles.

La loi de Dirichlet d'ordre K ≥ 2 de paramètres α1, ..., αK > 0 possède pour densité de probabilité :

f(x_1,\dots, x_{K}; \alpha_1,\dots, \alpha_K) = \frac{1}{\mathrm{B}(\alpha)} \prod_{i=1}^K x_i^{\alpha_i - 1}

pour tous les x1, ..., xK > 0 vérifiant x1 + ... + xK-1 < 1, où xK est une abréviation pour 1 – x1 – ... – xK–1. La densité est nulle en dehors de ce simplexe ouvert de dimension (K − 1).

La constante de normalisation est la fonction bêta multinomiale, qui s'exprime à l'aide de la fonction gamma :

\mathrm{B}(\alpha) = \frac{\prod_{i=1}^K \Gamma(\alpha_i)}{\Gamma\bigl(\sum_{i=1}^K \alpha_i\bigr)},\qquad\alpha=(\alpha_1,\dots,\alpha_K).

Propriétés[modifier | modifier le code]

Soit X = (X_1, \ldots, X_K)\sim\operatorname{Dir}(\alpha), signifiant que les K – 1 premières composantes possèdent la distribution précédente et que

X_K=1-X_1-\cdots-X_{K-1}.

Posons \textstyle\alpha_0 = \sum_{i=1}^K\alpha_i. Alors

 \mathrm{E}[X_i] = \frac{\alpha_i}{\alpha_0},

et

\mathrm{Var}[X_i] = \frac{\alpha_i (\alpha_0-\alpha_i)}{\alpha_0^2 (\alpha_0+1)},

En fait, les densités marginales sont des lois bêta :

X_i \sim \operatorname{Beta}(\alpha_i, \alpha_0 - \alpha_i).

Qui plus est,

\mathrm{Cov}[X_iX_j] = -\frac{\alpha_i \alpha_j}{\alpha_0^2 (\alpha_0+1)}.

Le mode de la distribution est le vecteur (x1, ..., xK) avec

 x_i = \frac{\alpha_i - 1}{\alpha_0 - K}, \quad \alpha_i > 1.


Agrégation[modifier | modifier le code]

Si X = (X_1, \ldots, X_K)\sim\operatorname{Dir}(\alpha_1,\ldots,\alpha_K),alors X' = (X_1, \ldots, X_i + X_j, \ldots, X_K)\sim\operatorname{Dir}(\alpha_1,\ldots,\alpha_i+\alpha_j,\ldots,\alpha_K). Cette propriété d'agrégation permet d'obtenir la distribution marginale de X_i mentionnée plus haut.

Distributions associées[modifier | modifier le code]

  • Si, pour \scriptstyle i\in\{1,2,\ldots,K\},
Y_i\sim \gamma(\alpha_i,1),, où γ désigne la distribution Gamma, indépendamment
alors
V=\sum_{i=1}^K Y_i\sim \gamma \left(\sum_{i=1}^K\alpha_i,1 \right),
et
(X_1,\ldots,X_K) = (Y_1/V,\ldots,Y_K/V)\sim \operatorname{Dir}(\alpha_1,\ldots,\alpha_K).
Bien que les Xi ne sont pas indépendants, ils peuvent néanmoins générer un échantillon de K variables aléatoires, distribuées selon une distribution Gamma. Malheureusement, puisque la somme V est perdue lors de la génération de X = (X1, ..., XK), il n'est pas possible de retrouver les variables Gamma initiales.

Génération de (pseudo-)nombres aléatoires (RNG)[modifier | modifier le code]

Une méthode pour obtenir un vecteur aléatoire x=(x_1, \ldots, x_K) à partir de la distribution Dirichlet de dimension K de paramètres (\alpha_1, \ldots, \alpha_K) est fournie par la remarque précédente. Tout d'abord, on tire K variables indépendantes y_1, \ldots, y_K selon des distributions Gamma, chacune avec la densité

 \textrm{Gamma}(\alpha_i, 1) = \frac{y_i^{\alpha_i-1} \; e^{-y_i}}{\Gamma (\alpha_i)}, \!

et on pose au final

x_i = y_i/\sum_{j=1}^K y_j. \!

Interprétations intuitives des paramètres[modifier | modifier le code]

Découpage d'une ficelle[modifier | modifier le code]

Une illustration de la distribution de Dirichlet apparaît lorsque l'on désire découper des ficelles (toutes de longueur initiale 1.0) en K pièces de différentes longueurs, où chaque pièce a, en moyenne, une longueur désignée mais cette longueur est autorisée à varier. Les valeurs α/α0 spécifient les longueurs moyennes des découpes résultant de la distribution. La variance (disparité autour de la moyenne) varie inversement avec α0.

Example of Dirichlet(1/2,1/3,1/6) distribution

Modèles d'urne et simulations du cas particulier des urnes de Pólya[modifier | modifier le code]

Considérons une urne contenant K différentes couleurs. Initialement, l'urne contient \alpha_1 boules de couleur 1, \alpha_2 boules de couleur 2, etc.. Procédons alors à N tirages dans l'urne suivant ce protocole: chaque boule tirée est replacée dans l'urne et on y ajoute une boule supplémentaire de même couleur. Lorsque N devient très grand, les proportions des boules de différentes couleurs sont distribuées selon \operatorname{Dir}(\alpha_1,\ldots,\alpha_K)[1].

Notons que chaque tirage modifie la probabilité d'obtenir une couleur donnée. Cette modification s'atténue d'ailleurs avec le nombre de tirages, puisque l'effet marginal d'ajout une boule supplémentaire diminue lorsque l'urne contient de plus en plus de boules.

Simulations de l'évolution de la proportion de boules blanches d'une urne de Pólya à deux couleurs[modifier | modifier le code]

Cas particulier d'une urne de Pólya à boules de deux couleurs : ce premier programme en R simule l'évolution de la proportion de boules blanches dans l'urne tout au long des n tirages.

# Une urne contient b boules blanches et r boules rouges
# On tire une boule au hasard et on la remet avec une boule supplémentaire
# de la même couleur, que l'on a pris dans la réserve.
# On note Xn le nombre de boules blanches au nième tirage.
#--polya3--Évolution de la proportion X/M de boules
#  blanches dans l'urne au cours de n tirages---
polya3 <- function(b = 8, r = 4, n = 100){
  urne <- rep(c("b", "r"), c(b, r))
  serieFreqb <- b / (b + r)
  for(k in 1:n){
    tiragek <- sample(urne, 1)
    if(tiragek == "r")     {
       urne <- c(urne, "r")} else {
       urne <- c(urne, "b")       }
    serieFreqb[k + 1] <- sum(urne == "b") / (b + r + k)
               }
# Affichages des résultats et graphiques
  plot(1:(n + 1), serieFreqb,
      pch = 21, col = "blue", bg = "blue", cex = .5,
      xlab = "Rang du tirage",
      ylab = "Proportion* simulée de boules blanches",
      main = paste("Évolution de la proportion* de boules blanches",
        "dans l'urne\n Simulée à partir d'une composition initiale de",
        b, "blanches et", r, "rouges"))
    abline(h = b / (b + r), col = "red")
}

Illustration des résultats de la fonction R : polya3() :

PolyaFreqRang.jpg


Simulations de la distribution de la proportion de boules blanches d'une urne de Pólya à deux couleurs[modifier | modifier le code]

Ce second programme en R simule la distribution de la proportion X/M de boules blanches de l'urne à la fin des n tirages.

#--polya2--Distribution simulée de la proportion X/M de boules
#  blanches dans l'urne au bout de n tirages---
polya2 <- function(b = 8, r = 5, n = 50, nbsim = 1000){
  seriePropX <- NULL ; classes <- seq(0, 1, .05)
  for(i in 1:nbsim){
    urne <- rep(c("b", "r"), c(b, r))
    for(k in 1:n){
      tiragek <- sample(urne, 1)
      if(tiragek == "r")     {
         urne <- c(urne, "r")} else {
         urne <- c(urne, "b")       }
                 }
    seriePropX[i] <- sum(urne == "b") / (b + r + n)
                   }
    moySerie <- mean(seriePropX)
    etSerie <- sqrt(var(seriePropX) * (nbsim - 1) / nbsim)
# Affichages des résultats et graphiques
  cat("Moyenne de la série des X/M  urne finale =", moySerie, "\n")
  cat("Ecart type de la série des X/M urne finale =", etSerie, "\n\n")
  par(mfrow = c(2, 1))
  hist(seriePropX, freq = FALSE, xlim = c(0, 1), breaks = classes,
      main = "Distribution simulée de la proportion finale de blanches",
      xlab = "Proportion finale de blanches dans l'urne")
    abline(v = b / (b + r), col = "red")
  boxplot(seriePropX, horizontal = TRUE, ylim = c(0, 1),
      main = paste("Effectifs initiaux : blanches =", b, "rouges =", r,
                   "\nnombre de tirages n=", n))
    abline(v = b / (b + r), col = "red")
}

Résultats numériques et graphiques de la fonction R : polya2() :

PolyaDistribFreq.jpg


Références[modifier | modifier le code]

  1. D. Blackwell and J. B. MacQueen 1973. Ferguson distributions via Pólya urn schemes. The Annals of Statistics, volume 1, number 2, pp353--355

Voir aussi[modifier | modifier le code]

Articles connexes[modifier | modifier le code]

Liens externes[modifier | modifier le code]

Non-Uniform Random Variate Generation, par Luc Devroye http://cg.scs.carleton.ca/~luc/rnbookindex.html