Revoir tout ce qui concerne les histogrammes et la programmation dans
Les histogrammes permettent de contrôler visuellement la justesse d'un
programme de simulation : il est indispensable d'en comprendre le
fonctionnement.
On indique, ici, quelques techniques générales de simulation de
variables aléatoires discrètes.
Dans le cas de lois connues on aura souvent avantage à utiliser de
méthodes particulières à ces lois !
La simulation de variables aléatoires est basée en Scilab
sur la fonction rand() (en fait, il en existe une autre
plus sophistiquée grand).
À chaque appel, cette fonction fournit un "tirage au hasard" :
un nombre, différent à chaque fois, compris entre 0 et 1
(strictement) suivant une loi de probabilité uniforme sur ]0,1[
et indépendant des autres tirages.
Pour simuler la réalisation ou non d'un événement de probabilité
p il suffit de se fixer un ensemble A de mesure de Lebesgue
p, contenu dans l'intervalle [0,1], et de considérer si
la valeur obtenue par l'appel de rand() tombe ou non dedans.
Bien entendu, le plus naturel est de prendre pour A un intervalle de
longueur p, par exemple [0,p] ou [1-p,1], ouvert ou
fermé, cela n'a aucune importance.
Si l'on attribue la valeur 1 à la réalisation de l'événement
et la valeur 0 à son contraire on dit que l'on a à faire à une
variable aléatoire de Bernoulli X. Ainsi,
p = 0.345;
X = double(rand() < p);
p.
La fonction double effectue la conversion %F --> 0 et
%T --> 1, mais elle n'est pas nécessaire lorsque l'expression booléenne
intervient dans une expression arithmétique comme dans
2*(rand(3,4) < p) - 1
3x4 aléatoire de +/-1.
Soit une variable X à valeurs dans l'intervalle d'entiers
E = 0:K dont la loi est donnée, par exemple par le vecteur
p = [ 0.4982, 0.2103, 0.1270, 0.0730, 0.0418, ..
0.0241, 0.0132, 0.0069, 0.0035, 0.0015, 0.0005 ];
K = length(p)-1sum(p) = 1 et
p(k+1) = P(X = k)
pour k de 0 à K
(les indices de vecteurs sous Scilab démarrent à 1, c'est
agaçant mais il faut en tenir compte !)
Pour simuler X on peut partitionner [0,1] en K+1
sous-ensembles Ak de mesures de Lebesgues
respectives p(k+1), on attribue, alors, à X la valeur
k lorsque rand() tombe dans Ak.
Bien entendu, le plus naturel (mais on peut procéder autrement, voir plus
bas) est de prendre la partition en intervalles consécutifs de longueurs
p(k) :
A0 = ]0,p(1)[, A1 = [p(1),p(1)+p(2)[,
A2 = [p(1)+p(2),p(1)+p(2)+p(3)[, etc.
On aimerait écrire quelque chose comme
U = rand(); // un tirage uniforme dans ]0,1[
if ( U < p(1) ) // 0 < U < p(1)
X = 0;
elseif ( U < p(1)+p(2) ) // p(2) <= U < p(1) + p(2)
X = 1;
elseif ...
N = 1600;
data = zeros(1,N);
for n = 1:N
U = rand();
k = 0;
F = p(1);
while ( U > F )
k = k + 1;
F = F + p(k+1);
end
data(n) = k;
end
s = -0.5:K+0.5;
histplot(s, data);
F (les valeurs de la fonction de
répartition) peuvent se calculer à l'avance une seule fois, et non 1600 fois,
en les stockant sous la forme d'un vecteur ordonné croissant.
De plus, on peut supprimer la boucle while
F = cumsum(p); // F = [p(1),p(1)+p(2),...,1]
N = 1600;
data = zeros(1,N);
for n = 1:N
U = rand();
data(n) = sum( U > F );
end
U > F est booléen, sum( U > F ) est le nombre
de fois où U > F.
On peut même supprimer la boucle for (sportif !)
F = cumsum(p);
N = 1600;
data = sum( rand(N,1)*ones(1,K+1) > ones(N,1)*F, 'c' );
F = cumsum(p);
N = 1600;
data = dsearch(rand(N,1), F);
dsearch fait une recherche dichotomique donc
très rapide dans le tableau ordonné F (il y a une petite
subtilité avec la valeur de retour 0).
--> help dsearch
Pour les variables aléatoires discrètes infinies on ne peut plus
utiliser les vecteurs p et F.
Il faudra définir la loi par une fonction et utiliser plutôt un code
analogue au premier exemple de programme (avec, cependant, l'avantage de
ne plus être enquiquiné par cette histoire de premier indice).
La semaine prochaine ?