Illustration du theorème central limit

Modélisation d'une loi discrète

La modélisation d'une loi discrète peut faire par une matrice à deux colonnes, la première représentant le support de la loi et la seconde la probabilité de tirage de l'évènement. Par exemple pour un dé à 6 faces :

d6 <- matrix( c(1:6,rep(1/6,6)), 6,2, dimnames=list(NULL, c('x','P')) )
d6
     x         P
[1,] 1 0.1666667
[2,] 2 0.1666667
[3,] 3 0.1666667
[4,] 4 0.1666667
[5,] 5 0.1666667
[6,] 6 0.1666667

Et le graphique correspondant :

plot(d6, pch=16, col='blue')

On peut maintenant suivant cette modélisation écrire une fonction qui calcule la loi d'une somme de deux lois de probabilités

somme <- function(loiA, loiB) {
    #Fonction générique mais qui gagnerait à être optimisée en utilisant une table de hash
    tailleMax <- length(loiA[,1]) * length(loiB[,1])
    support <- rep(0,tailleMax)
    probas <- rep(0,tailleMax)
    for (i in 1:length(loiA[,1])) {
        for (j in 1:length(loiB[,2])) {
            support[ (i-1)*length(loiB[,1]) + j ] = loiA[i,1] + loiB[j,1]
            probas[ (i-1)*length(loiB[,1]) + j ] = loiA[i,2] * loiB[j,2]
        }
    }
    index <- sort(support,index.return=TRUE)$ix
    r <- matrix(0,length(unique(support)),2,dimnames=list(NULL, c('x','P')))
    pos <- 0
    for (i in 1:length(support) ){
        if ((pos==0) || (value != support[index[i]]) ){
            value <- support[index[i]]
            pos <- pos + 1
            r[pos,1] <- value
        }
        r[pos,2] <- r[pos,2] + probas[index[i]]
    }
    return(r);
}
    
somme(d6,d6)  
      x          P
 [1,]  2 0.02777778
 [2,]  3 0.05555556
 [3,]  4 0.08333333
 [4,]  5 0.11111111
 [5,]  6 0.13888889
 [6,]  7 0.16666667
 [7,]  8 0.13888889
 [8,]  9 0.11111111
 [9,] 10 0.08333333
[10,] 11 0.05555556
[11,] 12 0.02777778  

Et le graphique de la loi correspondant à la somme de deux d6 :

plot(somme(d6,d6), pch=16, col='blue')

Au passage, l'on peut remarquer que la loi de la somme a en réalité été obtenue par un produit de convolution. Il existe des algorithmes permettant d'obtenir plus rapidement les résultats en passant dans l'espace de Fourrier.

Et maintenant on peut calculer la loi de la somme de 10 dés à 6 faces.

d <- d6
for (i in 1:9) {
    #on peut optimiser fortement en utilisant un arbre des sommes
    d <- somme(d,d6)
    }
plot(d, pch=16, col='blue')     

Je vous laisse supperposer avec une gaussienne de paramètres

m = 10*(7/2) = 35

et

sigma = sqrt(mean((d6[,1] - mean(d6[,1]))^2) * 10) = 5.400617

Pour mémoire la densité d'une loi normale (m,sigma) est donnée par

f(x) = 1 σ 2 π e - 1 2 x - m σ 2

Bref, pour une loi uniforme l'on converge très rapidement vers une loi normale, bien plus vite que ne l'aurait indiqué la borne de Markov.

Attention ce résultat suppose l'existence d'un moment d'ordre 2 (une variance) fini et est faux pour certaines lois classiques (Cauchy par exemple ou loi de student avec moins de 2 degrés de liberté).

On peut maintenant essayer en faisant varier la densité de probabilité initiale et le nombre de tirages additionnés et automatiser un peu le tracé du graphe.

afficheDensiteSomme <- function( loi , nbRep = 10 ){
    d <- loi
        for (i in 1:(nbRep-1)) {
        d <- somme(d,loi)
    }
    plot(d, pch=16, col='blue')
    m <- mean(loi[,1]) * nbRep
    s <- sqrt( mean((loi[,1]- mean(loi[,1]) )^2) * nbRep )
    x <- seq(min(d[,1]),max(d[,1]),by=0.1)
    lines(x,1/(s*sqrt(2*pi))*exp(-1/2*( (x-m)/s )^2) ,col='red')
}

loi <- d6
loi[6,1] <- 15 #la face 6 est remplacée par 15
afficheDensiteSomme(loi,5)

afficheDensiteSomme(loi,10)

afficheDensiteSomme(loi,30)

Design selector