vals=1:8
prob=c(0.1,0.2,0.05,0.05,0.05,0.15,0.35,0.05)
cdf=cumsum(prob)
#simulation de X:
x=runif(1)
which.max(x<cdf)
## [1] 2
#coût de génération = Espérance de X
ex=sum(vals*prob) #produit scalaire
# on optimise la génération en triant les probas à l'avance:
ps=sort(prob, decreasing=TRUE,index.return=TRUE)
p_decr=ps$x
cdf_decr=cumsum(p_decr)
ordre=ps$ix #il faut remettre les valeurs dans le bon "ordre"
# on genère avec:
x=runif(1)
ordre[which.max(x<cdf_decr)]
## [1] 6
#20 cases
#précalcul
nb_occurrences=prob*20
tableau=c()
for (i in 1:8){
tableau=c(tableau, rep(vals[i],nb_occurrences[i]))
}
#simulation
simtabulation=function(tableau){
sample(tableau,size=1,replace=TRUE)
}
simRejet=function(size,vals,prob){
x=1
y=2
while (y>prob[x]){ #reject couple (x,y)
x=sample(1:length(vals),size=size,replace=TRUE) #pick index
y=runif(1)*max(prob) #pick ordinate
}
vals[x] #if x values are different from 1:n
}
#on verifie la distribution empirique obtenue
N=1000
ech=c()
for(i in 1:N){
ech=c(ech,simRejet(1,vals,prob))
}
hist(ech)
Loi géométrique G(p) : nombre de tirages B(p) pour obtenir 1 succès
#loi geometrique
x=1
p=0.35
while (runif(1)<p){
x=x+1
}
x
## [1] 4
Loi binomiale Binom(n,p) : nombre de succès observés parmi n tirages B(p)
#loi binomiale (n,p)
n=50
sum(runif(n)<p)
## [1] 15
On peut toujours utiliser le rejet et l’inverse de la CDF.
Attention cependant pour la géométrique, le support est non borné, on ne peut pas utiliser l’algo de rejet tel quel.