Chapter 1Simulation de variables al´eatoiresR´ef´erences:[F] Fishman, A first course in Monte Carlo, chap 3.[B] Bouleau, Probabilit´es de l’ing´enieur, chap 4.[R] Rubinstein, Simulation and Monte Carlo Method, chap 3.Cadre: Scilab poss`ede une fonction rand() dont les appels successifs fournissent une suite devariables al´eatoires ind´ependantes et identiquement distribu´ees, de loi uniforme sur [0,1]. Nousne nous int´eresserons pas ici `a la conception d’une telle fonction.Probl`eme: Comment simuler une variable al´eatoire ou un vecteur al´eatoire suivant une loidonn´ee, diff´erente de la loi uniforme sur [0,1]?1.1 Premiers exemples de lois classiques1.1.1 Loi de Bernoulli de param`etre p: pδ +(1−p)δ1 0On propose deux fonctions. La premi`ere ne rend qu’une r´ealisation d’une variable al´eatoire dela loi de Bernoulli de parametre p, la seconde renvoie un k-´echantillon:function x=bernoulli1(p)// tirage suivant la loi de Bernoulli de parametre px=0; if rand()
Refe´rences: ´ [F] Fishman, A first course in Monte Carlo, chap 3. [B]Bouleau,Probabilite´sdel’inge´nieur,chap4. [R] Rubinstein, Simulation and Monte Carlo Method, chap 3.
Cadre:Scilabposse`deunefonction rand() dont les appels successifs fournissent une suite de variablesale´atoiresind´ependantesetidentiquementdistribu´ees,deloiuniformesur[0 1]. Nous nenousinteresseronspasicia`laconceptiond’unetellefonction. ´ Proble`me:Commentsimulerunevariableale´atoireouunvecteural´eatoiresuivantuneloi donn´ee,diff´erentedelaloiuniformesur[0 1]?
1.1 Premiers exemples de lois classiques 1.1.1LoideBernoullideparame`tre p : pδ 1 + (1 − p ) δ 0 Onproposedeuxfonctions.Lapremie`renerendqu’uner´ealisationd’unevariableale´atoirede la loi de Bernoulli de parametre p , la seconde renvoie un k -e´chantillon:
function x=bernoulli1(p) // tirage suivant la loi de Bernoulli de parametre p x=0; if rand()<p then; x=1; end;
function x=bernoulli2(k,p) // tirage d’un k-echantillon suivant la loi de Bernoulli de parametre p x=bool2s(rand(k,1)<p);
2
n 1.1.2Loibinomialedeparam`etres ( n p ) : X nk p k (1 − p ) n − k δ k k =0 Lemme 1.1.1 Si ( X i ) 1 ≤ i ≤ n sontdesvariablesale´atoiresind´ependantesetidentiquementdis-tribu´eesdeloideBernoullideparame`tre p , alors S = P in =1 X i suit une loi binomiale de param`etres ( n p ) .
function x=binomiale(n,p) // tirage suivant la loi binomiale de parametres (n,p) x=sum(bernoulli2(n,p)); 1.1.3 Loi uniforme sur { 0 1 n − 1 } : 1 n X − 1 δ k n k =0 function x=uniforme1(n) // tirage suivant la loi uniforme sur {0,...,n-1} x=floor(n*rand()); De´monstration: On note X lavariableal´eatoirerendueparcettefonction.Comme P (0 ≤ rand() < 1) = 1, on a P (0 ≤ n*rand() < n ) = 1 et P ( floor(n*rand()) ∈ [0 n − 1]) = 1 Donc X prend ses valeurs dans [0 n − 1]. Maintenant, soit k ∈ [0 n − 1]: P ( X = k ) = P ( floor(n*rand()) = k ) = P ( k ≤ n*rand() < k + 1) 1 1 = P ( kn ≤ rand() < k +) = n n
1.1.4 Loi uniforme sur [ a b ] : f ( x ) = b − 1 a 1 [ ab ] ( x ) Lemme 1.1.2 Si U suit la loi uniforme sur [0 1] , alors, si α > 0 et β ∈ R , αU + β suit la loi uniforme sur [ β β + α ] . De´monstration: Soit ϕ : R → R unefonctioncontinueborne´e:enfaisantle changement de variable x = αu + β , on a E ( ϕ ( αU + β )) = Z 01 ϕ ( αu + β ) du = Z ββ + α ϕ ( x ) α 1 dx ce qui signifie que X estunevariableale´atoirededensite´ f ( x ) = α 1 1 [ ββ + α ] , et on reconnaıˆtladensit´edelaloiuniformesur[ β β + α ]. 3
function x=uniforme2(a,b) // tirage suivant la loi uniforme sur [a,b] x=a+(b-a)*rand(); 1.1.5Loiexponentielledeparam`etre λ : f ( x ) = λ exp( − λx ) 1 R + ( x ) Lafonctionder´epartitiondelaloiexponentielleest F ( t ) = P ( X ≤ t ) = 1 − exp( − λt ) − Cette fonction est une bijection de ]0 + ∞ [ dans ]0 1[, d’inverse G ( u ) = λ 1ln(1 − u ) ♣ Exercice 1: Montrer que si U est de loi uniforme sur [0 1], alors G ( U ) suit la loi exponentielle deparame`tre λ .
function x=exponentielle(a) // tirage suivant la loi exponentielle de parametre a x=-log(rand())/a; ♣ Exercice 2: Pourquoia-t-onremplac´e1 − U par U dans l’algorithme? + ∞ 1.1.6Loig´eom´etriquedeparam`etre p : X (1 − p ) k − 1 pδ k k =1 M´ethode1:a`partirdujeudepileouface: Lemme 1.1.3 Si ( X i ) i ∈ N ∗ sontdesvariablesale´atoiresiiddeBernoullideparam`etre p , alors N = min { i : X i = 1 } suituneloige´ome´triquedeparametre p . `
function x=geometrique1(p) // tirage suivant la loi geometrique de parametre p x=1; while rand()>p, x=x+1; end; M´ethode2:`apartirdelaloiexponentielle: Lemme 1.1.4 Soit U suit la loi uniforme sur [0 1] , alors 1+ E lnl(n1 U − p ) suitlaloige´ome´trique deparame`tre p . D´emonstration: Notons X = 1 + E lnl(n1 U − p ) .Pard´efinitiondelapartieenti`ere, X prend ses valeurs dans N ∗ . Soit maintenant k ∈ N ∗ : P ( X = k ) = P 1 + E ln(l1n U − p ) = k = P E ln(l1n U − p ) = k − 1 = P k − 1 ≤ ln(l1n − Up ) < k = P ( k ln(1 − p ) < U ≤ ( k − 1) ln(1 − p )) = P (1 − p ) k < U ≤ (1 − p ) k − 1 = ≤ (1 − p ) k − 1 − (1 − p ) k = (1 − p ) k − 1 p
function x=poisson(a) // tirage suivant la loi de Poisson de parametre a test=exp(-a);x=0;prod=rand(); while (prod>=test), x=x+1; prod=prod*rand(); end; 1.1.8Loigaussiennecentr´eer´eduite:Me´thodedeBox-Mueller Rappelons la formule de changement de variable en plusieurs dimensions: The´ore`me1.1.6 Soit G un ouvert de R n , et g : G → R n une application de classe C 1 , injective et dont le jacobien det( J g ( x )) ne s’annule pas sur G . Soit f une application mesurable de g ( G ) dans R ,int´egrablesur g ( G ) . Alors Z f ( y ) dy = Z G f ( g ( x )) | det( J g ( x )) | dx g ( G ) Onend´eduitladensite´d’unvecteurimage: Proposition 1.1.7 Soit X unvecteurale´atoirea`valeurdans R n dedensit´e f X . Soit g : R n → R n une application de classe C 1 , injective et dont le jacobien det( J g ( x )) ne s’annule pas. Alors levecteural´eatoire Y = g ( X ) a pour densit´ e: f Y ( y ) = f X ( g − 1 ( y )) | det( J g − 1 ( y )) | D´ nstration: Soit ϕ : R n → R uneapplicationmesurable,positive,born´ee. emo Onappliqueleth´eor`emepr´ec´edentavec g − 1 : E ( ϕ ( Y )) = E ( ϕ ( g ( X ))) = Z R n ϕ ( g ( x )) f X ( x ) dx = Z R n ( y )) | det( J g − 1 ( y )) | ϕ ( y ) f X ( g − 1
cequiprouveler´esultat. En particulier, on obtient: Lemme 1.1.8 Soit R unevariableale´atoiredeloiexponentielledeparame`tre 1 2 et Θ une variableal´eatoiredeloiuniformesur [0 2 π ] ,suppose´esdeplusind´ependantes.Alors,sionpose X = R cos(Θ) et Y = R sin(Θ) ,lesvariablesal´eatoires X et Y sont iid de loi gaussienne centre´ere´duite.
function [x,y]=boxmueller // tirage de deux N(0,1) independantes r=sqrt(-2*log(rand())); t=2*%pi*rand(); x=r*cos(t); y=r*sin(t);