Si l'on compare à la distribution souhaitée:Simulons la variable aléatoire discrète X à valeurs dans {1,2,…,8} dont la distribution est fournie par le vecteur suivant:
valeurs_possibles=seq(1,8)
distri_probas=c(0.1,0.2,0.05,0.05,0.05,0.15,0.35,0.05)
On vérifie au passage que la somme des probas fait 1 :
sum(distri_probas)
## [1] 1
D'abord, on calcule la fonction de répartition, ou CDF (Cumulative Distribution Function), en calculant les sommes partielles des probabilités. Pour cela on écrit une fonction, qui sera réutilisable pour tout vecteur de probas.
repartition=function(p){
vec_rep=c()
vec_rep[1]=p[1];
for(i in 2:length(p))
vec_rep[i]=vec_rep[i-1]+p[i]
vec_rep
}
L'algo de simulation peut aussi s'écrire de façon générique par une fonction qui prend en argument le réel aléatoire uniforme sur [0,1] et le vecteur contenant la CDF.
simulation_cdf=function(x,v){
i=1
while(x>v[i]){
i=i+1
}
i
}
On aurait aussi pu faire les 2 en même temps (recherche et sommes partielles). Nous verrons cette alternative dans le cas de la loi géométrique, dont le support infini empêche de précalculer la fonction de répartition.
Essayons maintenant notre méthode :
set.seed(42) # au hasard...
simulation_cdf(runif(1),repartition(distri_probas))
## [1] 7
Bon, 1 échantillon, ça ne dit pas grand chose. Pour en faire un grand nombre:
N=100
fr=repartition(distri_probas)
echantillon_cdf=sapply(runif(N),simulation_cdf,fr) # applique la fonction simulation à tous les éléments de runif(N), avec comme second argument le vecteur fr
Combien a-t-on obtenu de 2?
sum(echantillon_cdf==2)
## [1] 17
Pas mal ! Vérifions graphiquement.
hist(echantillon_cdf, breaks=c(0,valeurs_possibles)) # on rajoute 0, sinon le calcul automatique des "boîtes" groupe les valeurs 1 et 2.
Si l'on compare à la distribution souhaitée:
barplot(distri_probas)
On remarque que les probabilités sont toutes multiples de 1/20. On peut donc remplir un tableau de 20 entiers choisis proportionnellement à ces probabilités, et piocher dedans uniformément.
nb_occ=(distri_probas)*20 # nb de répétitions de chaque valeur
tabul=c() #tableau des valeurs
for(i in 1:8){
tabul=c(tabul,rep(valeurs_possibles[i],nb_occ[i]))
}
echantillon_tabulation=sample(tabul,replace=TRUE,size=N)# prenons un echantillon de même taille qu'avant
hist(echantillon_tabulation, breaks=c(0,valeurs_possibles))
Rappel: la distribution binomiale \( Binom(n,p) \) vérifie : \( P(X=k)=C_n^k p^k (1-p)^{n-k} \). (On voit déjà venir les erreurs d'arrondi pour certaines valeurs de k…). Deux possibilités: on précalcule toutes les probas cumulées (comme pour l exercice 1) ou on les calcule à la volée.
binom_cdf=function(x,n=N,p=0.5){
seuil=p^N
k=0
while(x>=seuil){
k=k+1
seuil=seuil+choose(n,k)*p^(n-k)*(1-p)^k
}
k
}
stat_binom=sapply(runif(N*N),binom_cdf)
plot(stat_binom)
hist(stat_binom)
À comparer avec le générateur fourni par R:
echantillon_rbinom=rbinom(n=N*N,size=N,prob=0.5)
plot(echantillon_rbinom)
hist(echantillon_rbinom)
Pour la loi binomiale \( Binom(n,p) \), on peut aussi exploiter le fait que X représente le nombre de succès dans une succession de \( n \) tirages de Bernoulli \( \cal{B}(p) \)
binom_nb=function(n,p){
sum(runif(n)<p) # compte les coordonnées inférieures à p dans le vecteur de n nombres uniformes runif(n)
}
Ici le support est infini car la variable géométrique n'est pas bornée. On n'a pas le choix, il faut donc calculer la CDF à la volée.
geom_cdf=function(x,p=0.5){
seuil=p
k=1
while(x>seuil){
k=k+1
seuil=seuil+p^k
}
k
}
stat_geom=sapply(runif(N*N),geom_cdf)
plot(stat_geom)
hist(stat_geom)
Les puissances engendrent beaucoup d'erreurs d'approximation qui peuvent significativement biaiser l'échantillon. On peut obtenir une variable géométrique \( \cal G(p) \) comme le premier succès d'une succession (non bornée) de variables de Bernoulli de même paramètre \( p \).
geom_nb=function(p=0.5){
k=1
while(runif(1)>p){
k=k+1
}
k
}