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é optenue par un produit de convolution. Il existe des algorithmes permettant d'obtenir plus rapidement les résultat 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
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,30)
