Simularea variabilelor aleatoare

Generarea unor valori numerice n conformitate cu o anumita repartitie (simularea unei vari-
abile aleatoare) este o prelucrare frecvent utilizata atat n Statistica (de exemplu n tehnicile de
esantionare) cat si n modelarea stohastica si n implementarea algoritmilor aleatori (de tip Monte-
Carlo).
Exista doua modalitati principale de simulare a variabilelor aleatoare, si anume:
prin utilizarea unor dispozitive fizice cum ar fi zarul, ruleta si detectorii de radiatie nucleara;
prin implementarea pe calculator a unor algoritmi de generare a unor valori numerice.
Generatoarele bazate pe dispozitive fizice au cateva dezavantaje cum ar fi: rezultatele sunt
dificil de reprodus, iar functionarea dispozitivelor poate fi afectata de evenimente externe. Din
acest motiv, dupa aparitia si dezvoltarea calculatoarelor, interesul s-a orientat nspre metodele de
generare care se bazeaza pe capabilitatile aritmetice ale acestora.
Metodele de simulare a variabilelor aleatoare cu ajutorul calculatorului nu sunt perfecte ntrucat
calculatorul n sine este un dispozitiv determinist, iar majoritatea metodelor se bazeaza pe secvente
deterministe de calcul. Din acest motiv, valorile generate cu ajutorul calculatorului sunt n realitate
pseudo-aleatoare. In continuare vom omite calificativul pseudo doar pentru simplitate.
Din punct de vedere al modului de proiectare, generatoarele pot fi clasificate n:
generatoare de numere aleatoare uniform repartizate;
generatoare de numere aleatoare neuniform repartizate.
De cele mai multe ori generatoarele din a doua categorie se construiesc pe baza celor din prima
categorie.
In continuare vom prezenta algoritmi de generare a numerelor uniform repartizate si metode
generale pentru simularea variabilelor neuniform repartizate, discrete sau continui.

1 Simularea variabilelor aleatoare uniform repartizate

1.1 Cazul variabilelor discrete
O variabila aleatoare discreta, X, ce ia valori n multimea {1, . . . , n} este uniform repartizata daca
P (X = i) = 1/n pentru i = 1, n. Simularea unei astfel de variabile se bazeaza pe utilizarea unei
functii g : I k I, unde I este o submultime a numerelor ntregi care pot fi reprezentate n
calculator. Pentru a genera o secventa de numere se porneste de la cateva valori initiale x1 , . . . , xk
care formeaza samanta generatorului (seed) si se foloseste o relatie de recurenta:
xn = g(xn1 , . . . , xnk ), n>k (1)
Intrucat multimea valorilor numerice ce pot fi reprezentate n calculator este finita rezulta ca
sirul definit de relatia de recurenta (1) este periodic. Pentru ca generatorul implementat pe baza
acestei relatii sa fie acceptabil trebuie sa ndeplineasca cel putin doua conditii:

1
 Perioada lui sa fie mare n raport cu numarul de valori generate.

Valorile generate sa nu fie secvential corelate. Aceasta proprietate trebuie nteleasa n modul
urmator: considerand secvente de cate p valori succesiv generate se obtine o umplere uni-
forma a spatiului p-dimensional. In cazul prezentei corelatiei secventiale, aceste puncte din
spatiul p-dimensional se grupeaza ntr-un numar redus de hiperplane.
De remarcat ca importanta calitatii unui generator depinde de tipul aplicatiei unde este folosit
(ea este mai putin importanta n aplicatiile care utilizeaza numerele aleatoare doar pentru ilustrari
grafice de efect si devine foarte importanta de exemplu n aproximarea integralelor prin metoda
Monte-Carlo).
Satisfacerea conditiilor de mai sus poate fi asigurata printr-o buna alegere a functiei g. Cele mai
frecvent utilizate metode (pe care se bazeaza si generatoarele incluse n limbajele de programare)
sunt cele congruentiale. Acestea sunt caracterizate printr-o relatie de recurenta de forma:

xn = f (xn1 , . . . , xnk ) mod m

unde f : I k I este o functie iar k si m sunt valori care definesc generatorul.
Dintre metodele congruentiale, cele mai utilizate sunt cele pentru care f este o functie liniara:

xn = (a1 xn1 + . . . + ak xnk + c) mod m

cu a1 , . . . , ak , c si m valori ntregi care caracterizeaza generatorul. De regula ai , c {0, . . . , m 1}.
Evident relatia de recurenta se completeaza cu un set de valori initiale: x1 , . . ., xk .
Se observa ca aceste metode genereaza valori cuprinse ntre 0 si m 1, iar perioada poate fi cel
mult m. Din acest motiv m se alege cat mai mare, fiind de regula valoarea ntreaga maxima care
poate fi reprezentata n calculator. Mai dificila este nsa problema alegerii parametrilor a1 , . . . , ak ,
c. Pentru acestia nu exista retete riguroase, ci mai degraba exista valori particulare care au fost
determinate n urma a numeroase ncercari si testari ale eficientei.
Cel mai adesea se utilizeaza generatori de ordin 1, a caror forma este xn+1 = (axn + c)modm.
In cazul particular n care c = 0, generatorul se numeste congruential multiplicativ iar daca c 6= 0
e numit congruential mixt. Daca a si m sunt alese adecvat, atunci un generator multiplicativ este
la fel de bun ca unul mixt.
Cateva exemple de valori adecvate pentru parametrii unui generator de ordin 1 (k = 0) sunt:

1. a = 16807, c = 0, m = 231 1;

2. a = 1429 , c = 0, m = 231 1;

3. a = 24298, c = 99991, m = 199017;

4. a = 1993538837, c = 7261067085, m = 235 .

Secventa de valori generate depinde de valoarea initiala (x0 ). Aceasta este de regula stabilita
de catre utilizator prin proceduri specifice fiecarui limbaj de programare (de exemplu, procedura
randomize din Pascal sau functia srand din C). Daca se doreste ca la fiecare relansare a programului
sa fie generata alta secventa de valori este indicat ca x0 sa fie stabilita pornind de la ceasul intern
al calculatorului.
Daca m = 2 atunci se obtine un generator de biti conducand la o alta varianta de generare
a numerelor ntregi. Pentru o reprezentare pe 32 de biti, fiecare valoare ntreaga, xn , poate fi
exprimata ca:

2
 xn = b(1) (2) (32)
n bn . . . bn
(i)
iar fiecare bit bn se genereaza cu regula:
(i) (i)
b(i)
n = (bnp + bn(pq) )mod2

(i) (i)
n care p > q > 0, iar b0 , . . . , bp sunt valori binare care asigura initializarea generatorului.
Pentru ca generatorul sa aiba proprietati statistice bune, parametrii p si q se aleg astfel ncat
polinomul X p + X q + 1 sa fie prim n Z2 [X]. Un exemplu de alegere a parametrilor este p = 98 si
q = 27. In acest caz perioada generatorului este 2p 1.

1.2 Cazul variabilelor continui

Pentru a simula o variabila aleatoare uniform repartizata n intervalul [0, 1] este suficient sa se
genereze valori ntregi uniform repartizate pe multimea {0, . . . , m 1} si sa se mparta la m 1.
Daca se doreste simularea unei variabile uniform repartizate pe (0, 1) este suficient sa se genereze
valori n 1, . . . , m 1 si sa se mparta la m. De observat ca un generator mixt nu produce niciodata
valoarea 0.
Pe de alta parte, se considera ca un bun generator de valori uniform distribuite n intervalul
[0, 1] nu trebuie sa furnizeze niciodata valoarea 0 sau valoarea 1.
Pentru a simula o variabila, V , uniform repartizata n intervalul (a, b) este suficient sa se simuleze
o variabila, U , uniform repartizata n (0, 1), dupa care sa se faca transformarea V = (b a)U + a.
Majoritatea metodelor de simulare a variabilelor neuniforme se bazeaza pe un generator de valori
uniform repartizate n (0, 1). In cele ce urmeaza vom considera o metoda de generare a valorilor
uniform repartizate n (0, 1) ca fiind implementata de functia Random. Este usor de observat ca o
espresie de forma 1-Random simuleaza deasemenea o variabila uniform repartizata pe (0, 1).

1.3 Detalii de implementare

De remarcat ca unele dintre generatoarele liniar congruentiale de ordin 1 (exemplele 1, 2 si 4) nu pot
fi implementate n mod direct n limbaje de nivel nalt ntrucat valorile intermediare obtinute prin
calcule depasesc valorile maxime ce pot fi reprezentate ca ntregi pe 32 de biti. In aceste conditii,
majoritatea generatoarelor sunt implementate n limbaje de asamblare ce utilizeaza registrii de 64
de biti sau sunt utilizati algoritmi care produc doar valori intermediare ce nu depasesc 231 1.
In cazul simularii variabilelor continue, pentru a diminua efectul erorilor de rotunjire de la
operatia de mpartire este indicat sa se foloseasca reprezentarea n dubla precizie.
Intrucat generatorii liniar congruentiali produc valori n care cifrele de rang scazut (cele mai
putin semnificative) sunt mai putin aleatoare decat cele de rang nalt este indicat ca atunci cand se
doreste simularea unei variabile uniforme pe {1, . . . , s} sa se foloseasca catul impartirii si nu restul.
De exemplu, n C se va folosi o expresie de forma 1+(int)(s*rand()/(RAND MAX+1.0)) si nu una
de forma 1+(rand()%s).

2 Metoda inversarii functiei de repartitie

Fie X o variabila aleatoare (discreta sau continua) avand functia de repartitie FX : IR [0, 1]. Se
defineste inversa FX1 : [0, 1] IR n modul urmator:

3
 FX1 (u) = inf{x IR|FX (x) u}, u [0, 1].
Metoda inversarii functiei de repartitie se bazeaza pe urmatorul rezultat:
Fie F o functie de repartitie. Daca U este o variabila aleatoare uniform repartizata n
[0, 1], atunci functia de repartitie a variabilei aleatoare X = F 1 (U ) este F .
Intr-adevar, daca notam cu FX functia de repartitie a variabilei X = F 1 (U ) atunci are loc:
FX (x) = P ({X < x}) = P ({F 1 (U ) < x}) = P ({U < F (x)}) = P ({0 < U < F (x)}) = F (x)
deci functia de repartitie a lui X este chiar F .
Astfel forma generala a algoritmului de simulare a unei variabile avand functia de repartitie FX
este:
u:=Random
x:=InvF(u) /* InvF este inversa functiei F */
Return x
Pentru repartitii concrete se obtin variante particulare ale algorimilor, care difera ntre ele prin
modul de implementare a inversei functiei de repartitie.

2.1 Cazul repartitiilor discrete

Fie X o variabila aleatoare discreta care ia valori n multimea {x1 , x2 , . . . , xn }. Sa presupunem ca
x1 < x2 < . . . < xn , iar repartitia lui X este:
!
x1 . . . xi . . . xn
p1 . . . p i . . . p n
Functia de repartitie asociata va fi:


F1 = 0, x x1



F 2 = p1 , x1 < x x2



. . .
P
FX (x) = Fi = i1k=1 pk , xi1 < x xi



. . .

P

Fn = n1

k=1 pk , xn1 < x xn
Fn+1 = 1 x > xn
In acest caz valorile inversei functiei de repartitie pot fi calculate n modul urmator:
FX1 (u) = xi , daca FX (xi1 ) < u FX (xi ), pentru i = 1, n
cu x0 = si FX (x0 ) = 0.
Daca descriem functia de repartitie prin tabelul:
!
x0 = x1 . . . xi1 xi xi+1 . . . xn
F0 = 0 F1 . . . Fi1 Fi Fi+1 . . . Fn
atunci algoritmul de simulare consta n generarea unei valori, u, uniform repartizate n (0, 1) si n
gasirea indicelui i pentru care Fi1 < u Fi . Principala prelucrare a acestui algoritm este astfel
cautarea intervalului din lista valorilor functiei de repartitie care contine valoarea u.
Structura algoritmului este:

4
i:=1
u:=Random
While (u>Fi) and (i<n) Do i:=i+1
Return xi
Exemple.
1. Repartitia Bernoulli. Fie X : {0, 1} o variabila binara caracterizata de P ({X = 0}) = p
si P ({X = 1}) = q = 1 p. Functia de repartitie este:


0, x0
FX (x) = p, 0 < x 1

1, x > 1
In acest caz particular problema cautarii este mult simplificata si algoritmul devine:

u:=Random
If u<=p then x:=0 else x:=1
Return x

2. Repartitia binomiala. Fie X : {0, 1, . . . , n} o variabila avand functia de repartitie:

0, x0
FX (x) = P0 + . . . + Pk , k < x k + 1, k = 0, n 1

1, x>n
unde Pk = P (X = k) = Cnk pk (1 p)nk . Evident ca X poate fi simulata prin construirea
tabelului asociat functiei de repartitie si prin aplicarea metodei clasice de cautare. Pentru
valori mari ale lui n metoda nu este eficienta. O alta varianta este de a simula extragerile
cu revenire si de a numara de cate ori de produce evenimentul de probabilitate p. Pentru
aceasta se considera N ca fiind cel mai mic numar natural pentru care n1 = N p si n2 = N p
sunt si ele numere naturale. Se construieste un tabel cu n1 + n2 elemente avand urmatoarele
valori: t1 = . . . = tn1 = 1, tn1+1 = . . . = tn1+n2 = 0. Algoritmul pentru generarea unei valori
x este:

x:=0
For k:=1 to n do
{u:=Random
i:=INT(N*u)+1
x:=x+t[i]
}
Return x

3. Repartitia Poisson. Fie X : {0, 1, 2, . . .} o variabila aleatoare caracterizata prin

k
P (X = k) = e , cu > 0.
k!
Functia de repartitie este:
n1
X
FX (n) = P (X = k), n 1.
k=0

5
 In acest caz este ineficient sa se genereze un tabel al functiei de repartitie (ntrucat tabelul
este infinit si nu se cunoaste a priori cate elemente trebuie generate), fiind mai indicat ca evaluarea
acesteia sa se faca pe masura ce se efectueaza cautarea intervalului care contine valoarea u. Aceasta
varianta a algoritmului are forma:
u:=Random
i:=0; p:=exp(-lambda); f:=p
While (u>f) do {i:=i+1; p:=p*lambda/i; f:=f+p;}
Return i

2.2 Cazul repartitiilor continue

In cazul variabilelor continue exista doua situatii:
se cunoaste expresia analitica a inversei functiei de repartitie;
se foloseste un tabel cu valori ale functiei de repartitie.
In primul caz se foloseste forma generala a algoritmului, iar n al doilea se aplica algoritmul
specific repartitiilor discrete, nsa dupa ce s-a determinat intervalul (Fi1 , Fi ] care contine pe u se
calculeaza:

u Fi1
x = xi1 + (xi xi1 )
Fi Fi1
aceasta fiind valoarea care se returneaza si nu xi . In locul ultimei relatii, care asigura o interpolare
liniara a valorilor din tabelul functiei de repartitie, se poate aplica o formula corespunzatoare unui
alt tip de interpolare.
Exemple.
1. Repartitia exponentiala. Fie X : [0, ) o variabila aleatoare avand functia de repartitie:

FX (x) = 1 ex , > 0.
1
Inversa acestei functii este FX1 (u) = ln(1 u). Intrucat daca u este uniform repartizata

n (0, 1) atunci si 1 u este uniform repartizata n (0, 1), rezulta ca algoritmul de simulare
are forma:

u:=Random
Return -ln(u)/lambda

2. Repartitia Weibull. Fie X : [0, ) o variabila aleatoare avand functia de repartitie:

b
FX (x) = 1 eax , a, b > 0.
 1/b
1
Inversa lui FX este FX1 (u) = ln(1 u) astfel ca algoritmul de simulare este:
a
u:=Random
Return (-ln(u)/a)^(1/b)

Observatie. Metoda inversarii functiei de repartitie poate fi extinsa si n cazul variabilelor multi-
dimensionale. Pe o astfel de extindere se bazeaza metoda Box-Muller de simulare a variabilelor
normale. Aceasta metoda va fi prezentata n sectiunea dedicata variabilelor cu repartitie normala.

6
3 Metoda respingerii
Fie X o variabila aleatoare pentru care se cunoaste un algoritm de simulare si fie E un eveniment
de probabilitate nenula. Sa presupunem ca repartitia variabilei Y are proprietatea:

FY (y) = P ({Y < y}) = P ({X < y}|E). (2)

Astfel o valoare generata n conformitate cu repartitia lui X poate fi considerata o realizare a
variabilei Y daca este satisfacut evenimentul E. Realizarile lui X care nu satisfac evenimentul E
sunt respinse. Pornind de la aceasta idee se poate formula un algoritm general de simulare a lui Y :

Repeat
X:=GenerareX
Until E
Return X

Principalele probleme care trebuie rezolvate pentru a aplica metoda respingerii pentru simularea
unei variabile Y sunt:

gasirea variabilei X,

stabilirea evenimentului E,

astfel ncat sa fie satisfacuta relatia (2). In plus este de dorit ca n aplicarea algoritmului, numarul
mediu de respingeri ale valorilor generate sa nu fie prea mare. Numarul mediu de respingeri este
cu atat mai mic cu cat P (E) e mai apropiat de 1.
Un caz particular, destul de frecvent ntalnit, n care metoda respingerii poate fi aplicata cu
succes este cel al simularii repartitiilor uniforme multidimensionale pe domenii oarecare pornind de
la repartitii uniforme pe domenii paralelipipedice.
Presupunem ca Y , variabila care trebuie simulata, are repartitia uniforma pe domeniul D Rn .
Fie D0 = [a1 , b1 ] [a2 , b2 ] . . . [an , bn ] cel mai mic domeniu paralelipipedic care contine domeniul
D, iar X o variabila uniform repartizata pe D0 . Considerand evenimentul E = {X D}, variabila
Y poate fi simulata n modul urmator:

Repeat
/* generarea componentelor vectorului u din D */
For i:=1 to n Do ui=ai+(bi-ai)*Random
Until (u1, u2,..., un) in D

Se observa ca probabilitatea evenimentului E este vol(D)/vol(D0 ). Algoritmul va fi cu atat

mai eficient cu cat raportul ariilor este mai apropiat de 1. Pe de alta parte, complexitatea algo-
ritmului depinde de forma domeniului D ntrucat aceasta influenteaza complexitatea conditiei de
apartenenta.
Exemple.

1. Generarea unor puncte uniform repartizate n interiorul cercului unitate:

Repeat
x:=2*Random-1; y:=2*Random-1;
Until (x*x+y*y<1)

7
 2. Generarea unor puncte uniform repartizate n interiorul unui tor centrat n origine cu raza
mare egala cu 4 si cea mica egala cu 2:

Repeat
x:=-4+8*Random; y:=-4+8*Random; z:=-1+2*Random;
Until (z*z+Sqr(Sqrt(x*x+y*y)-3)<1)

Sa consideram problema mai generala a simularii variabilei Y avand densitatea fY pornind de

la variabila X cu densitatea fX . Presupunem ca exista c 1 astfel ncat

fY (x) cfX (x), x IR

Intuitiv metoda respingerii consta n acest caz n generarea de puncte uniform repartizate n
domeniul D0 delimitat de axa Ox si graficul functiei cfX si n acceptarea valorii abscisei doar daca
punctul apartine si domeniului D determinat de axa Ox si graficul functiei fY .
Astfel algoritmul de generare este:

Repeat
x:=GenerareX;
u:=Random;
Until c*f_X(x)*u<=f_Y(x);

Se observa ca evenimentul E este {cfX (X)U fY (X)} unde X este variabila aleatoare cu
densitatea fX iar U este o variabila uniform repartizata n (0, 1).
Aratam ca P (E) = 1/c. Intr-adevar:
Z
P (E) = P ({cfX (X)U fY (X)}) = fX (x)I(0,1) (u)dxdu =
{(x,u)|cfX (x)ufY (x)}

Z Z fY (x)/(cfX (x)) Z
fY (x) 1
= fX (x)dx I(0,1) (u)du = dx = .
{x|fX (x)>0} 0 {x|fX (x)>0} c c
Deci algoritmul este cu atat mai eficient cu cat c este mai apropiat de 1. Astfel este indicat
ca c sa se aleaga ca fiind cea mai mica valoare care verifica fY (x) cfX (x), de exemplu c =
maxx fY (x)/fX (x).
Pentru a justifica utilizarea
R
metodei trebuie sa aratam ca pentru orice multime boreliana, B,
are loc P ({X B}|E) = B fY (x)dx. Intr-adevar:
Z
P ({X B} E)
P ({X B}|E) = =c fX (x)I(0,1) (u)dxdu =
P (E) {xB}{cufX (x)fY (x)}

Z Z Z fY (x)/(cfX (x)) Z
=c fX (x)I(0,1) (u)dxdu = c fX (x) I(0,1) (u)du = fY (x)dx.
{xB}{ufY (x)/(cfX (x))} {xB} 0 B

Exemple.

1. Simularea repartitiei exponentiale trunchiata la (0, 1). Fie Y variabila avand densitatea de
repartitie: (
ex
1e 1 , daca x (0, 1)
fY (x) =
0 x 6 (0, 1)

8
 Fie X variabila cu densitatea de repartitie fX (x) = ex pentru x 0 si egala cu 0 n rest.
Atunci
fY (x) e
c = max = .
x0 fX (x) e1

Folosind pentru simularea variabilei X (cu repartitie exponentiala) algoritmul prezentat n

sectiunea anterioara, variabila Y poate fi simulata prin:

Repeat
x:=-ln(Random)
u:=Random
Until u<1/e

2. Simularea repartitiei seminormale. Fie Y o variabila cu densitatea de repartitie:

r
2 x2
fY (x) = e 2, daca x > 0

Consideram X o variabila cu repartitia exponentiala:

fX (x) = ex , daca x > 0.

Se observa ca: r
fY (x) 2e
c = max =
x fX (x)
deci algoritmul poate fi descris astfel:

Repeat
x:=-ln(Random)
u:=Random
Until ln(u)<=-(x-1)*(x-1)/2

3. Simularea abscisei unui punct uniform repartizat n discul unitate. Variabila de simulat, Y ,
are densitatea
2p
FY (x) = 1 x2 I[1,1] (x), xR

Considerand X ca fiind o variabila uniform repartizata n [1, 1] (fX (x) = 1/21[1,1] ) se
obtine c = 4/. Algoritmul poate fi scris n forma:

Repeat
x:=2*Random-1
u:=Random
Until (u*u<=1-x*x)

9
4 Metoda descompunerii
Se aplica atunci cand functia de repartitie a variabilei de simulat poate fi scrisa ca o combinatie
convexa a unor functii de repartitie asociate unor variabile pentru care se cunosc metode de simulare.
Fie F1 , F2 , . . . , Fn functii de repartitie si fie (p1 , p2 , . . . , pn ) o repartitie de probabilitate (0
P
pi 1, ni=1 pi = 1).
P
Consideram functia de repartitie F (x) = ni=1 pi Fi (x). Pentru a simula o variabila aleatoare n
conformitate cu repartitia data de F se aplica algoritmul:
Pas1. Se genereaza o valoare i {1, 2, . . . , n} n conformitate cu repartitia p (aplicand, de
exemplu, metoda inversarii functiei de repartitie).
Pas2. Se returneaza o valoarea generata prin metoda de simulare a variabilelor cu repartitia
corespunzatoare lui Fi
Observatie. Metoda se aplica n acelasi mod n cazul n care se cunosc functiile de densitate ale
variabilelor Xi .
Exemple.

1. Simularea repartitiei normale. Fie X o variabila cu repartitia normala standard. Atunci |X|
are repartitia seminormala. Tinand cont de faptul ca repartitia normala este simetrica putem
scrie:
1 1
FX (x) = F|X| (x) + F|X| (x),
2 2
astfel ca algoritmul de simulare a lui X este:

u:=Random
x:="generare variabila cu repartitia seminormala"
If u<0.5 then Return -x else Return x

2. Simularea variabilelor cu repartitie uniforma pe domenii decompozabile. Fie X o variabila

uniform repartizata pe D, un domeniu marginit. Presupunem ca D = D1 D2 . . . Dn
si Di Dj = pentru orice i 6= j. Pentru a aplica metoda descompunerii, se calculeaza
pi = ariaDi /ariaD si se aplica algoritmul:
Pas 1. Se genereaza i n conformitate cu (p1 , . . . , pn )
Pas 2. Se genereaza o valoare uniform repartizata n Di (utilizand de exemplu metoda
respingerii).

5 Simularea variabilelor aleatoare normal repartizate

5.1 Metoda bazata pe teorema limita centrala
Conform teoremei limita centrala, daca X1 , X2 , . . ., Xn sunt variabile aleatoare i.i.d. (independente
si identic repartizate) cu M (Xi ) = iar D2 (Xi ) = 2 atunci sirul de variabile aleatoare Zn =

n(X )/ tinde n distributie (cand n tinde la infinit) catre Z N (0, 1). X reprezinta (X1 +
. . . + Xn )/n.
Considerand X U( 0, 1) rezulta ca = 1/2, 2 = 1/12, deci

(X1 + . . . + Xn ) n/2
Zn = .
n/ 12

10
Se observa ca luand n = 12 expresia lui Zn devine foarte simpla (Z12 = X1 + . . . X12 6) astfel ca
un algoritm care simuleaza repartitia normala este:

z:=-6;
For i:=1 to 12 do
z:=z+Random
Return z

Observatie. Aceasta metoda este costisitoare din punct de vedere computattional ntrucat pentru
generarea unei valori z necesita 12 apeluri ale functiei Random.

5.2 Metoda Box-Muller

Metoda Box-Muller se bazeaza pe urmatorul rezultat teoretic:
Propozitie. Fie X N (0, 1) si Y N (0, 1) variabile independente si fie R si doua variabile
aleatoare ce au proprietatea ca X = Rcos() si Y = Rsin(). Atunci R si sunt independente, R
are repartitia Rayleigh (fR (r) = r exp(r2 /2) pentru r > 0 si 0 n caz contrar) iar are repartitie
uniforma pe [0, 2).
Din aceasta propozitie rezulta ca simuland o variabila R cu repartitia Rayleigh si una U(0,1)
putem genera doua valori n conformitate cu repartitia normala standard.
Simularea lui este evidenta, iar simularea lui R se poate face prin metoda inversarii functiei
de repartitie, caci:
q
FR (r) = 1 exp(r2 /2), r > 0 iar FR1 (s) = 2 ln(1 s).

Astfel algoritmul de generare poate fi scris n modul urmator:

u:=Random; v:=Random
r:=sqrt(-2*ln(u))
z1:=r*cos(2*pi*v)
z2:=r*sin(2*pi*v)
Return z1,z2

O varianta a acestui algoritm care evita folosirea functiilor trigonometrice si care se bazeaza pe
generarea de puncte uniform distribuite n interiorul cercului unitate este:

Repeat
u:=2*Random-1; v:=2*Random-1; s+u^2+v^2
Until 0<s<1
r:=sqrt(-2*ln(s)/s)
z1:=r*u; z2:=r*v
Return z1, z2

Observatie. La fiecare apel al oricareia dintre variantele algoritmului se obtin doua valori n con-
formitate cu repartitia normala standard.

11
6 Simularea variabilelor de tip 2 , Student si Fisher
Simularea acestor variabile se bazeaza pe legaturile ce exista ntre ele si repartitia normala stan-
dard. Presupunem ca avem la dispozitie o funtie normala care genereaza valori n conformitate cu
repartitia normala standard.
Pentru simularea unei variabile X 2 (n) se poate folosi algoritmul:

x:=0
For i:=1 to n do x:=x+normala(0,1)
Return x

Pentru simularea unei variabile X t(n) se poate folosi algoritmul:

z:=normala(0,1)
y:=0
For i:=1 to n do y:=y+normala(0,1)
x:=z*sqrt(n)/sqrt(y)
Return x

Pentru simularea unei variabile X F (n1 , n2 ) se poate folosi algoritmul:

y:=0
For i:=1 to n1 do y:=y+normala(0,1)
z:=0
For i:=1 to n2 do z:=z+normala(0,1)
x:=(y/n1)/(z/n2)
Return x

Observatie. Algoritmii folosesc din fiecare apel al al functiei normala ambele valori generate.

Referinte
[1] M.T. Boswell, S.D. Gore, G.P. Patil, C. Taillie, The Art of Computer Generation of Random
Variables, in Handbook of Statistics, vol. 9 (ed. C.R. Rao), Elsevier Science Publ., 1993.

[2] S.M. Ermakov, Metoda Monte Carlo si probleme nrudite, Ed. Tehnica, 1976.

[3] R.E. Crandall, Projects in Scientific Computation, Springer Verlag, 1994.

[4] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C. The Art
of Scientific Computing, Cambridge University Press, 1992.

[5] B. Ycart, Simulation des modeles markoviens, curs DESS dIngenierie Mathematique, Univ. J.
Fourier, Grenoble, 1997.

12