METODA MONTE-CARLO S ,I APLICAT ,II Conduc ator  stiint i c Conf. Dr. Tr^ mbit ,as,Radu Absolvent Costea Ionela Diana 2019 Cuprins Introducere 1… [618171]

UNIVERSITATEA BABES -BOLYAI CLUJ-NAPOCA
FACULTATEA DE MATEMATIC A S I INFORMATIC A
SPECIALIZAREA MATEMATIC A
LUCRARE DE LICENT  A
METODA MONTE-CARLO S ,I APLICAT ,II
Conduc ator  stiint ic
Conf. Dr. Tr^ mbit ,as,Radu
Absolvent: [anonimizat]
2019

Cuprins
Introducere 1
1 Preliminarii 2
1.1 Evenimete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Variabile Aleatoare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Valoarea medie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Dispersia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.5 Estimatori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.5.1 Estimatori nedeplasat ,i . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5.2 Estimatori corect ,i . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5.3 Estimatori absolut corect ,i . . . . . . . . . . . . . . . . . . . . . . 4
1.6 Inegalitatea lui Ceb^ as ,ev s ,i Legea numerelor mari . . . . . . . . . . . . . . 4
1.7 Variabile aleatoare de tip continuu . . . . . . . . . . . . . . . . . . . . . 5
2 Generarea numerelor aleatoare 9
2.1 Propriet at ,iile distribut ,iei de probabilitate . . . . . . . . . . . . . . . . . . 10
3 Metoda Monte Carlo de integrare 13
3.1 Integrarea Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Estimarea erorii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.3 Intervale de ^ ncredere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3.1 Intervale de ^ ncredere pentru media teoretic a . . . . . . . . . . . . 21
4 Tehnici de reducere a dispersiei 23
4.1 Metoda variabilei controlate . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Metoda select ,iei dup a important , a . . . . . . . . . . . . . . . . . . . . . . 29
4.3 Metoda select ,iei straticate . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.4 Metoda variabilei antitetice . . . . . . . . . . . . . . . . . . . . . . . . . 33
5 Aplicat ,ii 35
Bibliograe 48
0

Introducere
^In anul 1949 au fost puse primele baze teoretice ale acestei metode, s ,i anume metoda
Monte-Carlo, de c atre Ulam s ,i Metropolis. Dar ideea acestora a ap arut^ nca din anul 1777.
Atunci Buon a formulat faimoasa problem a a calculului probabilit at ,ii de intersect ,ie a
unui ac aruncat pe o suprafat , a plan a unde erau trasate drepte paralele s ,i echidistante.
Datorit a lui Barbier ^ n 1860, aceast a idee a permis calculul num arului prin arunc ari
succesive ale unui ac de lungime 1, pe o suprafat , a plan a cu drepte paralele la echidistanta
b. Metropolis, Fermi s ,i Ulam au folosit pentru aplicarea metodei, la cercetarea difuziei
neutronilor ^ n diferite materialele sionabile, liste de numere aleatoare care se publicau
pe atunci la Monte-Carlo. De aici vine s ,i cea mai cunoscut a denumire a metodei.
^In literatura de specialitate se ^ nt^ alnesc urm atoarele denumiri echivalente: metoda
^ ncerc arilor echivalente; metoda experiment arilor statistice; simularea numeric a; metoda
simularii indirecte; metoda numerelor aleatoare.
^Intr-o metod a Monte-Carlo, r aspunsul dorit este o cantitate ^ ntr-un model aleator s ,i
aceasta este estimat a de probe aleatoare ale modelului. Aplicat ,iile metodei Monte-Carlo
se pot folosi at^ at la estimarea integralelor c^ at s ,i la minimizarea funct ,iilor dicile sau
pentru a simula sisteme complexe.
Este o diferent , a important a^ ntre metoda Monte-Carlo s ,i metoda pseudo-Monte Carlo.
Prima estimeaz a cantit at ,i prin probe aleatoare iar cea de-a doua foloses ,te mostre care
sunt alese sistematic . ^Intr-un anumit sens, toate metodele practice de calcul sunt pseudo-
Monte Carlo, din moment ce generarea numerelor aleatoare implementat a ^ n mas ,ini nu
este cu adev arat la ^ nt^ amplare. Deci diferent ,a dintre cele dou a metode este put ,in neclar a.
Dar vom folosi termenul Monte Carlo pentru mostre care sunt generate folosind numere
pseudoaleatoare generate de un program computerizat. ^In cele ce urmeaz a ne vom axa
pe metoda Monte-Carlo.
Metodele Monte-Carlo sunt (cel put ,in ^ ntr-un anumit sens) metode ale ultimei sort ari.
Acestea sunt ^ n general destul de dicile s ,i se aplic a doar la problemele care sunt prea
greu de rezolvat cu metode deterministice.
1

1 Preliminarii
^In aceast a sect ,iune rezultatele s ,i not ,iunile teoretice sunt preluate din notele de curs
ale doamnei conf. dr. Natalia Ros ,ca [6].
1.1 Evenimete
Denit ,ie 1.1. Se numes ,te experiment orice proces sau act ,iune ale c arui rezultate nu
sunt cunoscute de dinainte cu certitudine, adic a sunt aleatoare.
Denit ,ie 1.2. Numim prob a ecare repetare a unui experiment. Iar orice rezultat posbil
al unui experiment se numes ,teeveniment .
Evenimetele elementare sunt evenimete care nu se pot realiza din cauza realiz arii
altor evenimente. Not am cu: !1;!2;:::evenimentele elementare. Fie
mult ,imea tuturor
evenimetelor elementare.

=f!1;!2;::;!ng
Atunci
se numes ,tespat ,iu de select ,ie.
Pentru oricare dou a evenimente As,iBdenim urmatoarele evenimente:
1.A[Bnumit reuniunea lui AcuB, adic a se realizeaz a cel put ,in unul dintre cele 2
evenimente.
2.A\Bnumit intersect ,ia luiAcuB, adic a se realizeaz a simultan ambele evenimete.
3.AnBnumitA-B,adic a se realizeaz a ev. As,i nu se realizeaz a ev. B.
4.A
Badic a se realizeaz a exact unul dintre evenimente.
Denit ,ie 1.3. Fie
o mult ,ime nit a. O familie nevid a de submult ,imiKP(
) se
numes ,te corp de p art ,i dac a:
(i)A2K=)A2K
(ii)A;B2K=)A[B2K
1.2 Variabile Aleatoare
Denit ,ie 1.4. Fie(
;K;P )c^ amp de probabilitate complet aditiv. O aplicat ,ieX:
!R
se numes ,tevariabil a aleatoare dac a:
8x2Ravem (X <x ) =f!2
jX(!)<xg2K
Clasicare :
2

1. Variabile aleatoare de tip discret: mult ,imea valorilor variabilei aleatoare Xeste cel
mult num arabil a.
2. Variabile aleatoare de tip continuu: mult ,imea valorilor variabilei aleatoare Xeste
un interval.
1.3 Valoarea medie
1. Cazul discret.
Denit ,ie 1.5. Fie variabila aleatoare Xde tip discret cu distribut ,iaX
xi
pi!
i2I.
Atunci spunem c a num arul real notat cu M(X)estevaloarea medie a lui Xs,i
este denit astfel:
M(X) =X
i2Ixi
pi
2. Cazul continuu.
Denit ,ie 1.6. Fie variabila aleatoare Xde tip continuu cu densitatea de probabil-
itatef. Atunci spunem c a num arul real notat cu M(X)estevaloarea medie a
luiXs,i este denit astfel:
M(X) =Z1
1x
f(x)dx
1.4 Dispersia
Denit ,ie 1.7. Fie variabila aleatoare Xcu valoarea medie m=M(X). Se numes ,te
dispersia variabilei aleatoare X num arul real notat cu D2(X), denit astfel:
D2(X) =M[(XM(X))2] =M[(Xm)2] (1.1)
Formula de calcul utilizat a ^ n practic a este:
D2(X) =M(X2)[M(X)]2(1.2)
1.5 Estimatori
Fie o select ,ie repetat a de volum n, X1;


;Xnvariabilele de selectie s ,ix1;


;xn
datele de select ,ie.
3

Denit ,ie 1.8. Numim estimator pentru parametrul necunoscut , statistica ^=^(X1;


;Xn)
care ia valori ^ n domeniul B R. Valoarea numeric a ^=^(x1;


;xn)se numes ,te
estimat ,ia lui
Denit ,ie 1.9. Estimatorul ^=^(X1;


;Xn)se numes ,teestimator consistent penru
parametrul dac a:
lim
n!1P(j^j<) = 1;8>0
1.5.1 Estimatori nedeplasat ,i
Denit ,ie 1.10. Estimatorul ^=^(X1;


;Xn)se numes ,teestimator nedeplasat
penru parametrul dac a:M(^) =
1.5.2 Estimatori corect ,i
Denit ,ie 1.11. Estimatorul ^=^(X1;


;Xn)se numes ,teestimator corect penru
parametrul dac a satisface urm atoarele condit ,ii:
a)lim
n!1M(^) =
b)lim
n!1D2(^) = 0
1.5.3 Estimatori absolut corect ,i
Denit ,ie 1.12. Estimatorul ^=^(X1;


;Xn)se numes ,teestimator absolut corect
penru parametrul dac a satisface urm atoarele condit ,ii:
a)M(^) =
b)lim
n!1D2(^) = 0
Propozit ia 1.13. Un estimator absolut corect este un estimator consistent.
Propozit ia 1.14. Un estimator corect este un estimator consistent.
1.6 Inegalitatea lui Ceb^ as ,ev s ,i Legea numerelor mari
Teorema 1.15. (Inegalitatea lui Ceb^ as ,ev)
Fie variabila aleatoare Xcu valoarea medie M(X)s,i dispersiaD2(X). Atunci avem
c a:
P(jXM((X)j<)1D2(X)
28>0 (1.3)
4

Adic aX2(M(X);M(X) +) cu o probabilitate 1D2(X)
2
Inegalitatea (1.3) se poate rescrie astfel:
P(jXM((X)j)D2(X)
28>0 (1.4)
Denit ,ie 1.16. Spunem c a s ,irul de variabile aleatoare (Xn)n1urmeaz a legea slab a a
numerelor mari dac a:
lim
n!1P 1
nnX
k=1Xk1
nnX
k=1M(Xk)<!
= 18>0 (1.5)
1.7 Variabile aleatoare de tip continuu
^In aceasta sect ,iune vom considera anumite tipuri de variabile aleatoare de tip
continuu.
Variabile aleatoare uniform distribuite
O variabil a aleatoare Xse spune c a urmeaz a Legea uniform a U(a;b) pe intervalul
(a;b),a<b , dac a densitatea de probabilitate este dat a de:
f(x) =8
<
:1
badac aa<x<b
0 ^ n rest
^In alte cuvinte, Xeste uniform distribuit ^ n intervalul ( a;b) dac a pune toat a masa ^ n
acel interval s ,i este probabil s a e "l^ ang a" orice punct din acel interval.
Media s ,i dispersia a unei variabile aleatoare uniform distribuit a ^ n intervalul ( a;b)
sunt obt ,inute astfel:
M(X) =1
baZb
axdx=b2a2
2(ba)=b+a
2
M(X2) =1
baZb
ax2dx=b3a3
3(ab)=a2+b2+ab
3
D2(X) =M(X2)[M(X)]2=1
12(ba)2
Prin urmare, valoarea as ,teptat a sau media este mijlocul intervalului ( a;b).
5

-3 -2 -1 0 1 2 300.050.10.150.20.250.30.350.4Figura 1: Funct ,ia de densitate
Funct ,ia de distribut ,ie a luiXeste dat a, pentru a<x<b , de:
F(x) =P(Xx) =Zx
a(ba)1dx=xa
ba
Variabile aleatoare normal distribuite
O variabil a aleatoare Xse spune c a urmeaz a Legea normal a N(;) cu media s,i
dispersia2dac a admite distribut ,ia:
f(x) =1p
2e(x)2
22;1<x<1
Densitatea normal a are forma unui clopot care este simetric fat , a de. (vezi gura 1).
Nu este greu de ar atat c a parametrii s,isunt egali cu media s ,i dispersia ^ n cazul
variabilelor aleatoare normal distribuite.Avem c a:
M(X) = D2(X) =2
Un lucru important despre variabilele aleatoare normal distribuite este c a, dac a X
este normal distribuit cu media s,i dispersia 2, atunci pentru orice constante as,ib,
aX+beste normal distribuit cu media a+bs,i dispersiaa22. Rezult a deci ca dac a X
6

−3 −2 −1 0 1 2 3−0.1−0.0500.050.10.150.20.250.30.350.4
−1.65
−1.268Figura 2:P(Z >za) =a
e normal distribuit media s,i dispersia2, atunci:
Z=X

este normal distribuit cu media 0 s ,i dispersia 1. O astfel de variabil a aleatoare Zse
spune c a are o distribut ,ie normal a standard. S a not am cu  funct ,ia de repartit ,ie a unei
variabile aleatoare normal a standard, atunci avem:
 =1p
2Zx
1ex2
2dx ;1<x<1
Acest rezultat este destul de folositor deoarece putem evalua toate probabilitat ,iile legate
deX^ n funct ,ie de . Spre exemplu, funct ,ia de repartit ,ie a luiXpoate  exprimat a
astfel:
F(x) =P(Xx) =P
X
x
!
=vP
Zx
!
= x

Valoarea lui ( x) poate  determinat a din tabel sau prin a scrie un program care s a
o aproximeze. Pentru a^ n intervalul (0 ;1),zas a e astfel ca:
P(Z >za) = 1(za) =a
Aceasta este, o normal a standard care va depas ,izacu probabilitatea a(vezi Figura 2).
7

Valoarea lui zapoate  obt ,inut a dintr-un tabel cu valorile lui . De exemplu, dac a
(1:64) = 0:95; (1:96) = 0:975; (2:33) = 0:99
observ am c a
z0:05= 1:64; z 0:025= 1:96; z 0:01= 2:33
Domeniul larg de aplicabilitate al variabilelor aleatoare normale rezult a din cea mai im-
portant a teorem a din teoria probabilitat ,ilor- Teorema Limit a Central a, care arm a c a
suma unui num ar mare de variabile aleatoare independente are aproximativ o distribut ,ie
normal a.
Teorema 1.17. Teorema Limit a Central a (Lindeberg-Levy)
Fie(ZN)N1un s ,ir de variabile aleatoare independente s ,i indentic distribuite (care
urmeaz a ecare aceeas ,i lege de probabilitate) cu:
M(ZN) =m D2(ZN) =2;8N1
unde
ZN=Zm
p
N;Z=1
NNX
k=1Zk
Atunci s ,irul(ZN)N1converge ^ n repartit ,ie c atre variabila alatoare Zcare urmeaz a legea
normal a N(0,1).
8

2 Generarea numerelor aleatoare
Vom genera numere aleatoare prin a lua mostre dintr-o colect ,ie de numere. Spre
exemplu, s a presupunem c a punem nc art ,i, numerotate de la 1 la n, ^ ntr-o cutie s ,i alegem
una la ^ nt^ amplare. Dup a ce ^ nregistr am num arul extras, punem cartea la loc. Am generat
un num ar aleator care este uniform distribuit printre valorile f1,2,…,ng, din moment ce
probabilitatea1
nde a extrage un anumit num ar este la fel pentru oricare alt num ar.
Dac a ad aug am kc art ,i ^ n cutie, toate numerotate cu 1, atunci distribut ,ia noastr a
devine neuniform a: probabilitatea de a alege 1 este acum(k+1)
(n+k)iar probabilitatea de a
alege orice num ar este1
(n+k).
Cele mai interesante sisteme aleatoare sunt acelea ale c aror rezultate sunt numere
dintr-un interval apart ,in^ and lui R2. Spre exemplu:
a) Dac a facem o roat a, prin punerea unui ac^ n centrul cercului. Desen am raza cercului.
^Inv^ artim acul, s ,i m asur am unghiul care ^ l formeaz a cu raza. Vom obt ,ine numere
aleatoare care sunt uniform distribuite ^ n intervalul [0 ;2).
b) Dac a, ^ n medie, o substant , a radioactiv a emite -particule la ecare secunde,
atunci timpul dintre dou a emisii succesive are distribut ,ia exponent ,ial a cu media 
(acesta este un caz special al distribut ,iei Gamma).
c) Distribut ,ia normal a, a c arui grac de probabilitate este sub forma unui clopot (vezi
gura 1), este un model bun in diferite situat ,ii:
-coecientul de inteligent , a IQ s-a construit astfel ^ nc^ at m asura s a e distribuit a
normal.
-multe caracteristici zice ale animalelor si plantelor (ex:^ nalt ,ie,greutate, etc) sunt
aproximativ normal distribuite.
-prezicerea probabilitat ,i ^ n jocurile de noroc a fost motivat ,ia lui DeMoivre(1667-
1754) pentru a deni distribut ,ia normal a.
Putem caracteriza un set de m asur atori aleatoare cu un set de numere (momente) derivate
din ele. Primele dou a momente, media s ,i dispersia, sunt cele mai importante. Valoria
medie a unui set de m asur atori fxig,i= 1;:::;n este
n1
nnX
i=1xi
9

s,i dispersia setului de probe m asoar a abatera probelor de la valoarea medie ca
n1
nnX
i=1(xi)2
Media s ,i dispersia setului de m asur atori sunt estimat ari ale mediei si dispersiei distibut ,iei
de la care au fost extrase. Dac a este cunoscut a media teoretica a distribut ,iei atunci dis-
persia distribut ,iei este estimat a de
n1
nnX
i=1(xi)2
2.1 Propriet at ,iile distribut ,iei de probabilitate
Sistemul particular din care sunt desenate monstrele aleatoare este caracterizat de
funct ,ia de densitate a lui nenegativ a f(x). Domeniul
al funct ,iei este mult ,imea tu-
turor valorilor posibile care ar putea  obt ,inute dac a lu am monstre aleatoare ale funct ,iei.
Intervalul funct ,iei este o submult ,ime a numerelor nenegative. Valoarea
Z
f(x)dx
este probabilitatea ca o mostr a din aceast a distribut ,ie este ^ n mult ,ime denit de
(dac a domeniul este discret atunci integrala se ^ nlocuies ,te cu o sum a). Datorit a acestor
propriet at ,i, putem observa c aZ
f(x)dx= 1
(^In cazul discret, suma probabilitat ,ilor este 1.)
Spre exemplu, pentru cele ncart ,ii din cutia de mai sus, domeniul este
= f1;2;3;:::;ng
s,if(x) =1
npentrux= 1;2;:::;n . Pentru roata noastr a, domeniul este
= [0 ;2) s,i
f(x) =1
2.
Mediaa distribut ,iei(cunoscut a ca s ,i valoarea as ,teptat a) este suma tuturor
valorilor posibile, luate dup a probabilitatea lor. ^In alte cuvinte, este valoarea as ,teptat a
care am putea s a o calcul am f ac^ and media valorilor obt ,inute pentru media mostrelor
dac a am lua un num ar mare de mult ,imi foarte mari de probe. Similar, dispersia2
distribut ,ieieste estimat a de dispersiile mostrelor. Media si dispersia distribut ,iei sunt
denite prin :
=Z

xf(x)dx
10

2=Z

(x)2f(x)dx
Din nou, pentru cazul discret, ^ nlocuim integrala cu o sum a. Deci pentru cele nc art ,i
din cutie avem:
=nX
i=1i1
n=n+ 1
2
2=nX
i=1
in+ 1
2!2
1
n=(n+ 1)(n1)
12
Vom da funct ,ile de distribut ,ie pentru 3 distribut ,ii continue ment ,ionate ^ n exemple:
Distribut ,ia uniform a pentru intervalul [0 ;m] are funct ,ia de densitate:
f(x) =1
m
Media si dispersia sunt :
=mZ
0x
mdx=m
2
2=mZ
01
m
xm
22
dx=m2
12
Distribut ,ia exponent ,ial a cu parametrul are
f(x) =1
ex

pentru x2[0;1). Media este iar dispersia este 2=2.
Distribut ,ia normal a cu parametrul s,iare
f(x) =1p
22e(x)2
22
pentrux2(1;+1). Media este iar dispersia este 2.
Exemplu: [6] S a se calculeze media M(X) s,iM(X2) pentru Legea Binomial a:
X0
B@k
Ck
npkqnk1
CA
cup+q= 1.
11

Rezolvare: Utiliz am formula de recurent , a pentru combin ari:
Ck
n=n
kCk1
n1=n(n1)
k(k1)Ck2
n2
S,i formula binomului lui Newton :
nX
k=0Ck
npkqnk= (p+q|{z}
=1)n= 1
M(X) =nX
k=0k
pk=nX
k=0k
Ck
n
pk
qnk=nX
k=1k
Ck
n|{z}
=n
Ck1
n1
pk
qnk=n
nX
k=1Ck1
n1
pk
qnk
=n
pnX
k=1Ck1
n1
pk1
qnkNOT :i=k1=n
pn1X
i=0Ci
n1
pi
qn1i=
=n
p
(p+q|{z}
=1)n1=n
p
=>M (X) =n
p
M(X2) =nX
k=0k2
pk=nX
k=0k2
Ck
n
pk
qnk=nX
k=1k2
Ck
n
pk
qnk=
=nX
k=1[k2k+k]
Ck
n
pk
qnk=nX
k=1k(k1)
Ck
n
pk
qnk+nX
k=1k
Ck
n
pk
qnk
|{z}
=M(X)=
=nX
k=2k
(k1)
Ck
n|{z}
n(n1)Ck2
n2
pk
qnk+M(X)|{z}
=n
p=n
(n1)nX
k=2Ck2
n2
pk
qnk+np=
=n
(n1)
p2nX
k=2Ck2
n2
pk2
qnk+npNOT :i=k2=n
(n1)
p2n2X
i=0Ci
n2
pi
qni2
|{z}
= (p+q)n2+np=
=n
(n1)
p2
(p+q|{z}
=1)n2+np=n
(n1)
p2+np
=>M (X2) =n
(n1)
p2+np
12

3 Metoda Monte Carlo de integrare
Multe metode estimeaz a excelent valoare integralei
I=Z

f(x)dx
atunci c^ and feste o funct ,ie neted a de o singur a variabil a xs,i
este un interval nit [ a;b].
Spre exemplu, putem ^ mp art ,i intervalul ^ n subintervale mici, construim o aproximare
polinomial a pentru f^ n ecare subinterval prin evaluarea funct ,iei la mai multe puncte
din acel subinterval, aproxim am integrala peste subinterval cu integrala polinomial a s ,i
^ nsum am valorile aproximative. Aceasta este ideea din spatele metodei Newton-Cotes.
Pentru integralele multidimensionale situat ,ia este mai dicil a. S a presupunem c a
x2R10iar
este hipercubul unitate [0 ;1]


 [0;1].Atunci, ca s ,i exemplu, un polinom
de gradul 2 ^ n ecare variabil a va avea termeni de forma:
x[]
1x[]
2x[]
3x[]
4x[]
5x[]
6x[]
7x[]
8x[]
9x[]
10
unde num arul din ecare chenar este 0,1 sau 2.
Deci polinomul are 310= 59;049 coecient ,i, adic a ne-ar trebui 59,049 de valori ale
funct ,iei pentru a rezolva aceast a problem a. Este o metod a destul de complicat a.
O alternativ a ar  s a aproxim am f(x) printr-o funct ,ie separabil a:
f(x)f1(x1)f2(x2)


f10(x10)
Desigur, nu putem face acest lucru pentru orice funct ,ie. Dar dac a acest lucru este posibil
atunci putem aproxima integrala prin:
IZ1
0f1(x1)dx1


Z1
0f10(x10)dx10
Aceast a metod a nu merge de ecare dat a. De aceea ne-ar trebui alt a opt ,iune s ,i anume
13

Metoda Monte Carlo pentru integrare.
Algoritm 1: Estimarea volumului prin metoda Monte Carlo
Fie Rdo mult ,ime care cont ,ine
s ,i are volumul Jcunoscut. Gener am n
punctez(i)2 care sunt uniform s ,i aleator distribuite. Fie ^ nnum arul acestor punctelor
care de asemenea se a
 a ^ n
s ,i estim am volumul lui
prin vnJunde
vn=^n
n:
Algoritm 2: Integrarea Monte Carlo
Gener amnpunctez(i)2
care sunt uniform s ,i aleator distribuite. Atunci valoarea
medie a lui f^ n
este aproximat a de:
n=1
nnX
i=1f(z(i))
iar valoare integralei este:Z

f(x)dx=InZ

dx
Sunt dou a metode simple de a folosi metoda Monte Carlo pentru integrare. Pentru
prima, s a ne g andim la integrare ca s ,i cum am calcula volumul unui solid
, s ,i ^ ncorpor am
acest solid ^ ntr-unul mai mare pentru care calculul volumului este us ,or. Gener am o
succesiune de puncte aleatoare din s ,i determin am fract ,iavde puncte care apart ,in de
asemenea lui
. Putem estima volumul lui
prin v^ nmult ,it cu volumul lui . As ,a
obt ,inem Algoritmul 1.
Ca s ,i exemplu trivial consider am un sfert dintr-un cerc
pe care ^ l ^ ncorpor am ^ ntr-un
p atrat de 0,5 pentru a putea calcula aria.
14

Figura 3: Estimarea ariei 0.19636 a unui sfert de cerc folosind metoda Monte Carlo cu
20 de evaluari de functii
Alt punct de vedere ne ofer a un algoritm mai bun, unul prin care obt ,inem mai multe
informat ,ii de la valorile funct ,iei. Ideea este simpl a: conform teoremei valorii medii a
integralei, integrala
I=Z

f(x)dx
este egal a cu valoarea medie a lui fpe
^ nmult ,it cu volumul lui
, adic aR

dx. Acesta
ne ofer a a dou a metod a, Algoritmul 2.
C^ at de bun a este aproximarea lui Inde la Algoritmul 2? Valoarea as ,teptat a a lui
Ineste valoarea adev arat a a integralei. Deasemenea, pentru un nmare, estim arile au
aproximativ o distribut ,ie normal a cu dispersia2
n, unde
2=Z

(f(x)n)2dx
Iar dispersia este o constant a care depinde de variat ,ia ^ nf^ n jurul valorii lui medii,
dar nu de dimensiunea da domeniului de integrare
.
15

3.1 Integrarea Monte Carlo
Problemele care se rezolv a, ^ n general, prin metoda Monte-Carlo sunt acelea ale
estim arii valorii medii a unei variabile aleatoare. Cazul c^ and momentul de ordinul al doilea
al variabilei aleatoare este nit este cel mai us ,or: calcul am N realiz ari independente ale
variabilei aleatoare iar valoarea medie a acestei variabile o estim am cu ajutorul mediei
aritmetice a acestor realiz ari. Deoarece dintre tot ,i estimatorii nedeplasat ,i ai valorii medii,
media aritmetic a are cea mai mic a dispersie acest procedeu se poate studia ^ n am anunt.
Fief:Rd!Ro funct ,ie,Xun vector aleator iar f(X) o variabli a aleatoare cu
densitatea de probabilitate :Rd!R. Atunci avem c a:
M(f(X)) =Z
Rdf(x)(x)dxNOT:=I (3.1)
undeM(X) este valoarea medie a variabilei aleatoare f(X). Este necesar a integrabili-
tatea p atratului lui f(x). Av^ and ^ n vedere c a feste liniar a pe spat ,iul considerat deducem
c aR
Rdf(x)2(x)dx< +1.
Ne propunem s a estim am integrala (3.1). Pentru aceasta vom genera Nnumere
aleatoare s ,i independente xi,i= 1;


;N. FieX1;


;XNvariabile aleatoare inde-
pendente s ,i identic distribuite pe intervalul [0 ;1]dcu funct ,ia de densitatea obis ,nuit a.
Atunci avem c a:
J=1
NNX
k=1f(Xk)
iar valoarea medie este:
J=1
NNX
k=1f(xk)
unde Jeste media variabilelor aleatoare f(X1);


;f(XN):
Propozit ia 3.1. ([5]) Estimatorul Jare urm atoarele propriet at ,i:
1.M(J) =I(estimator nedeplasat al lui I)
2.D2(J) =D2[f(X)]
N
16

3.lim
N!1D2(J) = 0
4.P
lim
N!1J=I
= 1
Observ am c a estimatorul Jsatisface condit ,iile 1 s ,i 3, deci este un estimator absolut
corect pentru I. Putem concluziona c a formula de integrare a metodei Monte Carlo este:
I=Z
Rdf(x)(x)dxJ=1
NNX
k=1f(xk) (3.2)
Observat ,ie 3.2. Dac a(x) = 0 ,x =2DiarDRsatunci integrala din formula (3.1)
devine:
I=Z
Df(x)(x)dx:
Unde variabilele aleatoare X1;X2;


;XN, sunt independente s ,i identic distribuite ^ n D.
3.2 Estimarea erorii
Dac a valoarea estimat a nu este egal a cu valoarea as ,teptat a nu spunem c a este
o gres ,eal a, ci doar un rezultat al variat ,iei numerelor aleatoare. Variabilitatea datelor
aleatoare se observ a^ n eroarea experimental a. Ca^ n orice problem a statistic a de estimare,
valoarea estimat a ar trebui s a e ^ nsot ,it a de valoarea estimat a de c atre dispersia sa.
Estimarea dispersiei a unui estimator de interes este, de obicei, doar dispersia es ,antionului
de valori calculate ale estimatorului de interes.
Una dintre metodele folosite pentru estimarea erorii^ n metoda Monte Carlo se bazeaz a
pe inegalitatea lui Ceb^ as ,ev. Reamintim formula (1.3):
P(jXM((X)j<)1D2(X)
28>0
Not am cu2(f) abaterea standard a variabilei aleatoare f(X) s,i rescriem formula
(1.1):
2(f) =D2[f(X)] =M(f(X)2)[M(f(X))]2;(cazuldiscret ) (3.3)
17

2(f) =D2[f(X)] =Z
Rsf(x)2(x)dx"Z
Rsf(x)(x)dx#2
;(cazulcontinuu ) (3.4)
Propozit ia 3.3. ([3]) ^In formula de integrare Monte Carlo (3.2) avem urm atoarea esti-
mare a erorii:
P 1
NNX
k=1f(Xk)Z
Rsf(x)(x)dx<(f)p
N
!
1
;8
2(0;1) (3.5)
Demonstrat ie: Aplic am estimatorului Jinegalitatea lui Ceb^ as ,ev:
P(jJM(J)j<)1D2(J)
2;8>0 (3.6)
Din propozit ,ia 3.1 avem c a:
M(J) =I D2(J) =2(f)
N
^Inlocuim ^ n formula (3.6) s ,i obt ,inem:
P 1
NNX
k=1f(Xk)I<!
12(f)
N
2;8>0
Alegem acum =(f)p
N

=)P 1
NNX
k=1f(Xk)Z
Rsf(x)(x)dx<(f)p
N
!
1
;8
2(0;1)
O alt a modalitatea pentru estimarea erorii ^ n metoda Monte Carlo se bazeaz a pe
Teorema limit a central a (1.17). Folosind notat ,iile din Teorem a putem arma c a ZN=Z,
adic a variabila aleatoare ZNpoate  considerat a c a urmeaz a legea normala N(0,1), c^ and
N!1:
Teorema 3.4. ([3]) Estimarea erorii ^ n formula de intergare Monte Carlo (3.2), c^ and
18

N!1 este dat a de:
P 1
NNX
k=1f(Xk)Z
Rsf(x)(x)dx(f)p
N!
=2()1;8>0; (3.7)
unde este funct ,ia lui Laplace, funct ,ia de distribut ,ie a legii normale N(0,1), adic a:
() =1p
2Z
1ex2
2dx:
Demonstrat ie: Aplic am Teorema Limit a Central a s ,irului de variabile aleatoare ( f(XN))N1,
iar acestea sunt independente s ,i identic distribuite. Atunci avem c a:
M[(f(XN)] =I; D2[f(XN)] =2(f); N1:
ZN=1
NNX
k=1f(Xk)I
(f)p
N=JI
(f)p
N
Atunci s ,irul de variabile aleatoare ( ZN)N1urmeaz a legea normal a N(0,1), cum N!
1, conform Teoremei Limite Centrale. Fix am >0 s,i atunci putem scrie:
P(jZNj)=2()1:
^InlocuimZNs,i obt ,inem:
P
jJIj(f)p
N!
=2()1:
Observat ,ie 3.5. ^In cazul ^ n care abaterea standard (f)este necunoscut a, putem estima
parametrul 2(f)prin estimatorul nedeplasat:
2(f) =1
N1NX
k=1(f(Xk)J)2: (3.8)
19

3.3 Intervale de ^ ncredere
Fie caracteristica X cu legea de probabilitate f(x;), unde parametrul 2BR
este necunoscut. Se consider a o select ,ie repetat a de volum niarX1;


;XNvariabilele
de select ,ie s ,ix1;


;xNdatele de select ,ie. Fie num arul 2(0;1), numit probabilitate
de risc iar 1, numit probabilitate de ^ ncredere .
Denit ,ie 3.6. [6] Se numes ,teinterval de ^ ncredere pentru parametru , intervalul
aleator (1;2) = ( 1(X1;


;XN);2(X1;


;XN)), unde statisticile 1s,i2satisfac
condit ,ia:
P(2(1;2)) = 1
Intervalul numeric (1;2) = ( 1(x1;


;xN);2(x1;


;xN))se numes ,tevaloarea in-
tervalului de ^ ncredere pentru parametrul necunoscut 
Propozit ia 3.7. [6] Fie caracteristica Xce urmeaz a legea normal a N(m;)s,i o select ,ie
repetat a de volum n. FieX1;


;Xnvariabilele de select ,ie. Atunci statictica
T=Xm
p
N(3.9)
urmeaz a legea Student T(n1), cun1grade de libertate.
Demonstrat ie: Not am cu
U=Xm
p
NV=1
2nX
k=1(XkX)2;
unde U urmeaz a legea normal a N(0;1) iar V legea 2cun1 grade de libertate.
Ar at am c a
T=Uq
V
n1:
20

Avem c a:
Uq
V
n1=Xm
p
N
pn1
1
pPn
k=1(XkX)2=Xmp1
n1Pn
k=1(XkX)2
pn=Xm
pn=T
3.3.1 Intervale de ^ ncredere pentru media teoretic a
Not ,iunile teoretice din aceast a sect ,iune sunt preluate din [6]
Fie caracteristica X ce urmeaz a legea normal a de parametrii ms,i, undem2R
(media teoretic a) este necunoscut s ,i>0 este deasemenea necunoscut. Fie statistica
T=Xm
p
N
Din propozit ,ia (3.4) s ,tim c a aceasta urmeaz a legea Student T(n1).
Determin am intervalul numeric ( t1;t2) astfel ^ nc^ at:
P(T2(t1;t2)) =Fn1(t2)Fn1(t1) = 1
undeFn(x) este funct ,ia de repartit ,ie a legii Student cu ngrade de libertate denit a astfel:
Fn(x) =(n+1
2))pn(n
2)Zx
1
1 +t2
n!n+1
2
dt; x2R
Aleg^ andt2=t1
2s,it1=t2, avemFn1(t1
2) = 1
2, iarP( m1<m< m2) = 1.
Atunci intervalul de ^ ncredere pentru media teoretic a are extremit at ,ile:
m1=Xt1
2
pn
m2=X+t1
2
pn
Putem concluziona c a un interval (1 )% de ^ ncredere pentru I=R
Rdf(x)(x)dx
21

este:
Jt1
2
(f)p
N;J+t1
2
(f)p
N!
(3.10)
22

4 Tehnici de reducere a dispersiei
Din formulele (3.5) s ,i (3.7) putem remarca faptul c a estimarea erorii ^ n metoda
Monte Carlo, mai exact limitele erorii, depind dep
Ns,i(f).
Putem cres ,te ecacitatea metodei prin diverse modalit at ,i legate de dispersie. Aceste
tehnici de reducere a dispersiei dateaz a din 1953, de la ^ nceputul calculelor digitale. Dar
sunt la fel de importante s ,i ^ n zilele noastre.
Reducerea dispersiei este de fapt, c autarea unei alternative s ,i a unor estimatori mai
precis ,i ai unei cantit at ,i date.
Posibiltatea de a reduce dispersia este ceea ce separ a metoda Monte Carlo de simularea
direct a. Tehnicile simple de reducere a dispersiei sunt, ^ n general, remarcabil de eciente
s,i us ,or de implementat. ^In anumite apliicat ,ii, precum ^ n Chimia Cuantic a, reus ,esc ceea
ce ar  imposibil ^ n caz contrar.
Ideea principal a este de a folosi informat ,ii suplimentare legate de funct ,ia de integrat
s,i de domeniul pe care aceasta se integreaz a, ^ n ordine s a reduc a efectul de es ,antionare
aleatoare. Acesta este unul dintre principiile fundamentale ale proiect arii statistice.
Dintre toate metodele de reducere a dispersiei, care pot  folosite s ,i ^ n combinat ,ie, am-
intim: metoda select ,iei dup a important , a, metoda variabilei controlate, metoda variabilei
antitetice s ,i metoda select ,iei straticate.
^In continuare vom prezenta c^ ateva metode de reducere a dispersiei.
4.1 Metoda variabilei controlate
^In aceast a metod a vom folosi covariant ,a. Orice variabil a care este corelat a cu
variabila de interes are valoarea potent ,ial a pentru reducerea dispersiei estimatorului.
O asemenea variabil a este folositoare dac a este us ,or de generat s ,i dac a are propiet at ,i
cunoscute sau pot  calculate us ,or.
Fie X o variabil a aleatoare iar metoda Monte Carlo introduce estimarea mediei =
M(X). Presupunem c a Y este o variabil a aleatoare cu media cunoscut a M(Y), atunci
23

e variabila aleatoare:
~X=X+b[YM(Y)]8bconst:
este un estimator nedeplasat pentru . Pentru a determina cea mai bun a valoare a
luib, folosim covariant ,a:
D2[X+b(YM(Y))] =D2(X+bY) =D2(X) +b2D2(Y) + 2bCov (X;Y )
Pentru a reduce dispersia, valoarea optim a a lui b este
bNOT:=Cov(X;Y )
D2(Y):
Iar pentru aceast a valoare dispersia estimatorului este:
D2(X+b(YM(Y))) =D2(X)[Cov(X;Y )]2
D2(Y)(4.1)
Cantitatea Y este numit a variabil a controlat a pentru simularea estimatorului X. Pen-
tru a vedea dac a funct ,ioneaz a, observ am c a beste negativ (pozitiv) c^ and X s ,i Y sunt
corelate pozitv (negativ). S a presupunem c a X s ,i Y sunt corelate pozitiv, adic a X este
mare c^ and Y este mare s ,i vice versa. Dac a rezultatul unei simul ari produce o valoare
mare (mic a) a lui Y (este indicat de faptul c a Y este mai mare (mic) dec^ at media lui
cunoscut a) atunci este probabil adev arat c a X este de asemenea mai mare (mic) dec^ at
media lui, s,i am dori s a corect am acest lucru prin a reduce (cres ,te) valoarea estima-
torului X, iar acesta este facut din moment ce beste negativ(pozitiv). Analog cazului
c^ and X s ,i Y sunt corelate negativ.
^Imp art ,ind ecuat ,ia (4.1) laD2(X) obt ,inem:
D2(x+b(YM(Y))
D2(X)= 1Corr2(X;Y )
24

unde coecientul de corelat ,ie dintre X s ,i Y este:
Corr (X;Y ) =Cov(X;Y )p
D2(X)D2(Y)
^In generalCov(X;Y ) s,iD2(Y) sunt necunoscute s ,i trebuie estimate din datele simu-
late. Dac a consider am nsimul ari cu datele obt ,inuteXi;Yi;i= 1;


;n, atunci folosind
estimatorii
dCov(X;Y ) =nX
i=1(XiX)(YiY)
n1(4.2)
s,i
bD2(Y) =nX
i=1(YiY)2
n1(4.3)
putem aproxima bprin^b, unde
^b=Pn
i=1(XiX)(YiY)Pn
i=1(YiY)2:
Dispersia estimatorului controlat
D2(X+b(YM(Y))) =1
n
D2(X)Cov2(X;Y )
D2(Y)!
poate  estimat a folosind estimatorul covariant ,eiCov(X;Y ) ^ mpreun a cu estimatorii
dispersieiD2(X) s,iD2(Y).
Cazul concret Rescriem (3.1) ^ ntr-o form a echivalent a:
I=Z
Rsg(x)(x)dx+Z
Rs[f(x)g(x)](x)dx: (4.4)
Alegem funct ,iag:Rs!Rastfel ca prima integral a din formula (4.4) s a e integrabil a.
Atunci aceasta va  us ,or de calculat analitic. Cea de a doua integral a poate  estimat a
cu ajutorul metodei Monte Carlo. Dac a funct ,ia g copiaz a comportamentul funct ,iei f,
astfel ^ nc^ at s a absoarb a aproape toat a dispersia funct ,iei f, atunci vom obt ,ine o reducere
25

a dispersiei.
Estimarea integralei init ,iale se va reduce atunci la estimarea integralei:
I1=Z
Rs[f(x)g(x)](x)dx:
Folosind metoda Monte Carlo aplicat a anterior, estimatorul integralei I1este:
IA=1
NNX
k=1(f(XK)g(Xk)]
Not am cu
Y=f(X)g(X);
unde X este un vector aleator cu funct ,ia de densitate .
Comparat ,ia dispersilor lui IAs,iJdin formula (3.2) se reduce la comparat ,ia dispersiei
lui Y cu dispersia lui f(X). Acestea satisfac urm atoarea relat ,ie:
D2(Y) =D2(f(X)) +D2(g(X))2Cov(f(X);g(X))
=)D2(Y)<D2(f(X))()D2(g(X))<2Cov(f(X);g(X)):
Vericarea condit ,iei de mai sus o vom face experimental, folosind punctele xk; k=
1


;Ns,i formulele (4.2) s ,i (4.3):
D2(g) =1
N1NX
i=1"
g(Xi)1
NNX
k=1g(Xk)#2
Cov(f;g) =1
N1NX
i=1"
f(Xi)1
NNX
k=1f(Xk)#"
g(Xi)1
NNX
k=1g(Xk)#
Dac a condit ,ia:D2(g)<2Cov(f;g) este satisf acut a ^ nseamn a c a funct ,ia g a fost sucient
de bine aleas a.
26

Exemplu: ([9]) (Blackjack) Acest foc de c art ,i este, ^ n general, jucat cu dealer-ul care
amestec a mai multe pachete de c art ,i, pune de o parte c art ,ile folosite s ,i ^ n nal reamestec a
c art ,ile c^ and num arul de c art ,i r amase este mai mic dec^ at o anumit a limit a. S a pre-
supunem c a o rund a nou a ^ ncepe de ecare dat a c^ and dealer-ul reamestec a c art ,ile s ,i
suntem interesat ,i s a folosim simul ari pentru a estima M(X), c^ as ,tigurile unui juc ator pe
o rund a, unde ne asum am faptul c a juc atorul aplic a o anumit a strategie x a, de exemplu
s a numere c art ,ile deja folosite ^ n aceea rund a s ,i mizeaz a diferite sume ^ n funct ,ie de acea
"num ar atoare". Presupunem c a jocul const a ^ ntr-un juc ator ^ mpotriva dealer-ului.
Procesul aleator rezult a din amestecarea c art ,ilor de c atre dealer. Dac a acesta foloses ,te
kpachete a c^ ate 52 de c art ,i, putem genera amestecarea cu ajutorul unei permut ari
aleatoare de numere de la 1 p^ an a la 52 k. Not am cu I1;


;I52kaceste permut ari. Dac a
x am
uj=Ijmod 13 + 1
s,i l as am
vj=min(uj;10)
atuncivj;j= 1


;52kreprezint a valorile succesive ale c art ,ilor amestecate cu 1 ind
notat asul.
Fie N num arul de m^ aini jucate ^ ntr-o tur a, s ,iBjsuma pariat a pe m^ ana j. Pentru
a reduce dispersia, putem folosi o variabil a controlat a care s a e mare c^ and juc atorul
det ,ine mai multe m^ aini bune dec^ at deler-ul s ,i este mic a ^ n cazul contrar. Din moment ce
a avea 19 sau mai mult este bine, denim urm atoarele:
Wj= 1 dac a a 2-a carte a juc atorului pus a pe pariul jface s a e 19 iar Wj= 0 ^ n
cazul contrar.
Zj= 1 dac a a 2-a carte a dealer-ului pus a pe pariul jface s a e 19 iar Zj= 0 ^ n cazul
contrar.
Este evident c a Wjs,iZjau aceeas ,i distribut ,ie, atunciM(WjZj) = 0 s ,i nu este
27

dicil de ar atat c a:
M"NX
j=1Bj(WjZj)#
= 0:
Putem folosi ca s ,i variabil a controlat aPN
j=1Bj(WjZj). Bine^ nt ,eles, nu e sigur ca
19 este cea mai bun a valoare. Dar sunt c^ ateva lucr ari preliminare care indic a faptul c a
19 funct ,ioneaz a cel mai bine, s ,i a rezultat ^ n reducerea dispersiei de 15% sau chiar mai
mult, depinde de strategia folosit a de juc ator.
O reducere mai mare a dispersiei ar trebui s a e dac a folosim dou a variabile controlate.
Una ind denit a ca mai sus, cu except ,ia c aWjs,iZjsunt sigur 1 dac a m ana este 19 sau
20. Iar cea de a doua variabil a este similar a, doar ca de aceast a dat a indicatorii sunt 1
c^ and este blackjack ^ n m^ an a.
C^ and folosim mai multe variabile controlate calculele pot  f acute cu ajutorul unui
program de calculator pentru regresii liniare multiple.
X=a+kX
i=1ciYi+e
undeeeste o variabil a aleatoare cu media 0 s ,i dispersi2. L as ambb
is a e estimarea celui
mai bunbi; i= 1;


;katunci avem c a:
bb
i=bci; i = 1


;k
undebci; i = 1;


;ksunt cele mai mici p atrate ale regresiei estimate de ci. Valoarea
estimatei controlate poate  obt ,inut a astfel:
X+kX
i=1bb
i(YiM(i)) =ba+kX
i=1bciM(i)
Aceasta este estimat ,ia controlat a care este de fapt estimat ,ia regresiei liniare multiple
avaluate ^ n punctul ( M(1);


;M(k)).
Dispersia estimat ,iei controlate poate  obt ,inut a prin ^ mp art ,irea regresiei lui 2cu
28

num arul de simul ari rulate.
4.2 Metoda select ,iei dup a important , a
Aceast a metod a se numes ,te as ,a deoarece se bazeaz a pe as ,a numitele funct ,ii im-
portante, distribut ,iile esent ,iale, care t ,in locul distribut ,iilor originale. De asemenea, se
bazeaz a pe o reprezentare alternativ a a formulei (3.1). O bun a alegere a noii distribut ,ii
ne v-a ajuta la reducerea dispersiei, lucru pe care ni-l dorim.
Fie atunciDRss,ifunct ,ia de densitate strict pozitiv a.
Aceast a metod a, select ,ia important a, ne sugereaz a s a alegem o nou a funct ,ie de den-
sitateh:D!Rcare este strict pozitiv a pe acest domeniu. Atunci integrala noastr a I
se va rescrie astfel:
I=Z
Df(x)(x)dx=Z
Df(x)(x)
h(x)h(x)dx: (4.5)
Folosind metoda Monte Carlo, denim estimatorul ISIal integralei I din formula (4.5)
astfel:
ISI=1
NNX
k=1f(Yk)(Yk)
h(Yk)(4.6)
undeYk,k= 1;


;Nsunt variabile cu funct ,ia de densitate h.
Not am acum cu Yvectorul aleator care are, de asemenea, funct ,ia de densitate h.
Atunci avem c a:
D2(ISI) =1
ND2"
f(Y)(Y)
h(Y)#
=1
NZ
Dh
f(x)(x)
h(x)Ii2
h(x)dx
Comparat ,ia dispersilor lui ISIs,iJdin formula (3.2) se reduce la comparat ,ia dispersiei
luif(Y)(Y)
h(Y)cu dispersia lui f(X), unde X este un vector aleator cu funct ,ia de densitate
.
29

S,tim c a:
D2"
f(Y)(Y)
h(Y)#
=M"
f2(Y)2(Y)
h2(Y)#
"
M
f(Y)(Y)
h(Y)!#2
Obiectivul acestei metode este, deci, s a alegem h^ n as ,a fel ^ nc^ at dispersia aceasta s a e
minimalizat a. Deoarece
"
M
f(Y)(Y)
h(Y)!#2
= Z
Df(x)(x)dx!2
Alegerea implic a doar primul termen din expresia dispersiei. Folosind inegalitatea lui
Jensen, obt ,inem o limit a inferioar a pentru acel termen:
M"
f2(Y)2(Y)
h2(Y)#
M"
jf(Y)jj(Y)j
h(Y)#2
= Z
Djf(x)(x)jdx!2
Este evident atins a c^ and
h(x) =jf(x)(x)jR
Djf(x)(x)jdx
4.3 Metoda select ,iei straticate
^In aceast a metod a, anumite proport ,ii ale domeniului de integrare sunt luate din
regiuni specicate ale spat ,iului. Unul dintre obiective este s a ne asigur am c a toate regiu-
nile sunt acoperite. Un alt obictiv este s a reducem dispersia general a a estimatorului prin
select ,ii mai bune unde funct ,ia "dispare", adic a valorile f(xi) pot varia foarte us ,or.
Select ,ia straticat a se efectueaz a, ^ n general, prin formarea unor subregiuni distincte
cu funct ,ii diferite s ,i importante ^ n ecare dintre acestea. Se bezeaz a deci pe alegerea unor
puncte din regiunile ^ n care funct ,ia nu dispare.
Consider am DRss,i X un vector aleator cu valori ^ n acest domeniu D. Fiefunct ,ia
30

de densitate a acestui vector. Prin aceast a metod a vom estima valoare integralei:
I=Z
Df(x)(x)dx (4.7)
Vom reprezenta domeniul de integrare sub forma unei reuniuni nite de subdomenii
disjuncteDk; k = 1;


;m. Avem atunci c a:
I=Z
Df(x)(x)dx=mX
k=1Z
Dkf(x)(x)dx: (4.8)
S a consider am acum restrict ,iile funct ,iei din formula (4.8) la subdomeniile introduse
mai susDk,8k= 1;


;m, denite astfel:
fk:Dk!R; fk(x) =f(x);8x2Dk
Fiecare din integralele din membrul drept al ecuat ,iei (4.8) se estimeaz a cu ajutorul
metodei Monte Carlo. Atunci vom deni:
gk=Z
Dk(x)dx;  k(x) =8
><
>:(x)
gkx2Dk
0x =2Dk
=)Ik=Z
Dkf(x)k(x)dx:
Cu aceste not ,iuni introduse formula (4.8) devine:
I=mX
k=1gkIk
Pentru reducerea dispersiei, metoda select ,iei straticate presupune estimarea ec arei
integraleIk. Fix am atunci n1;


;nmnumere ^ ntregi astfel ^ nc^ at n1+


+nm=n.
Gener amnk;k= 1


;mvariabileXk;1;


;Xk;nk, cu aceas ,i distribut ,ie s ,i funct ,ia de
31

densitatek. Denim urm atoarea variabil a aleatoare astfel :
Sk=1
nknkX
i=1f(Xk;i)
Observ am c a:
M(Sk) =Z
Dkf(x)k(x)dx=Ik
Estimatorul select ,iei straticate I , notat cu ISSeste astfel:
ISS=mX
k=1gkSk: (4.9)
Deoarece
M(ISS) =mX
k=1gkM(Sk) =I:
adic a satisface condit ,ia (1.5.1), este un estimator nedeplasat. Iar pentru c a
D2(ISS) =mX
k=1g2
kD2(Sk) =mX
k=1c2
k1
nk"Z
Dkf2(x)k(x)dxI2
k#
putem concluziona pe baza condit ,iilor (1.5.3) c a ISSeste un estimator absolut corect
pentru integrala I.
Teorema 4.1. ([7]) Fie date n s ,iDk;k= 1


;msubdomenii. Atunci dispersia estima-
torului select ,iei straticate, ISSeste minimalizat de alegerea ca nks a e de forma:
nk=gkkPm
l=1cll
n; k = 1


;m
unde2
k=D2(f(Xk;1)):
Teorema 4.2. ([7]) Dac ank=ngkpentruk= 1;


;m, atunci avem urm atoarea relat ,ie
^ ntre dispersia estimatorului de select ,ie straticat a ISSs,i dispersia estimatorului metodei
32

standard Monte Carlo J:
D2(ISS) =D2(J)1
nmX
k=1gk
Ik
gkI!2
:
Observat ,ie 4.3. ([7]) ^In practic a, este dicil de aplicat teroema (4.1) deoarece 1;


;m
sunt ^ n general necunoscute. Iar teorema (4.2) poate  destul de greu de aplicat, ^ n
cazurile ^ n care este dicil de generat vectorul aleator X, condit ,ionat deX2Dk.
4.4 Metoda variabilei antitetice
Ne dorim s a calcul am integrala I din formula (3.1) s ,i avem la dispozit ,ie o funct ,ie
g:D!R;sucient de apropiat a de f(x) astfel ca
I=Z
Df(x)(x)dx=Z
Dg(x)(x)dx (4.10)
Obt ,inem urm atorul estimator al integralei (4.10) folosind metoda Monte Carlo:
IAN=1
NNX
k=11
2[f(Xk) +g(Xk)]
Iar estimat ,ia este:
IAN=1
NNX
k=11
2[f(xk) +g(xk)]
undexk;k= 1


;Nsunt puncte aleatoare av^ and funct ,ia de densitate , iarXk;k=
1


;Nsunt variabilele corespunz atoare.
Putem observa c a avem nevoie de 2 Nevalu ari ale funct ,iei. Din aceast a cauz a vom
avea nevoie de estimatorul care foloses ,te 2Nevalu ari ale funct ,ieif, obt ,inut prin metoda
Monte Carlo astfel:
J=1
2N2NX
k=1f(Xk): (4.11)
Comparat ,ia dispersilor lui IANs,iJdin formula (4.11) se reduce la comparat ,ia dis-
33

persiei estimatorului 1cu dispersia estimatorului 2, denite astfel:
1=1
2[f(X1) +g(X1)]
2=1
2[f(X2) +f(X2)]
Observ am c a D2(2) =1
2D2[f(X)], unde X este un vector aleator. Atunci putem
scrie:
D2(1) =1
4D2[f(X1) +g(X1)]
=1
4D2[f(X)] +1
4D2[g(X)] +1
2Cov(f(X);g(X))
=1
2D2(2) +1
4"
D2[g(X)] + 2Cov(f(X);g(X))#
Dac aD2[g(X)] + 2Cov(f(X);g(X))<0, atunci avem
D2(1)<1
2D2(2)<D2(2):
C^ andf(X) s,ig(X) sunt negativ corelate, obt ,inem o alegere buna pentru g deoarece
D2[g(X)]>0:
34

5 Aplicat ,ii
(1)S a se calculeze ^ n R folosind metoda Monte Carlo urm atoarea integral a:
Z1
0 Z1
x Z2
xycos(xy)ezdz!
dy!
dx
Figura 4: Domeniul de integrare
f <- function(x,y,z) cos(x*y)*exp(z);
intMC1 <- function(N){
inte <- 0;
i=0;
while (i < N) {
x <- runif(1,0,1); y <-runif(1,0,1); z <- runif(1,0,2);
if (x <=y & x*y <= z) {
i <- i+1;
inte = inte +f(x,y,z);} }
inte <- inte/N*2;}
int <- c(0,0,0); intA <- c(0,0,0);
for (k in 1:3) {
int[k]=intMC1(1e4);
35

Figura 5: aplicat ,ia 1
intA[k]=intMC1(1e6);}
int
intA
z<-intMC1(1e7)
z
Rezultatul afis ,at:
> int
[1] 6.586625 6.608632 6.598312
> intA
[1] 6.589457 6.595819 6.591859
> z<-intMC1(1e7)
> z
[1] 6.592951
>
36

(2)S a se calculeze ^ n R folosind metoda Monte Carlo volumul si masa unui obiect dat
s,tiind c axyz1,5x5,5y5 iar5z5:
Z Z Z
(x;y;z )dxdydz;
unde(x;y;z ) =e0:5z
Figura 6: Domeniul de integrare
rho <- function(z) exp(0.5*z);
N<-10000000;
volumeOfBox<-1e3;
vol<-0; mass<-0;
volsq<-0; masssq<-0;
i <- 0;
while (i < N) {
x<-runif(1,-5,5); y<-runif(1,-5,5); z<-runif(1,-5,5);
if (x*y*z <= 1) { i <- i+1;
37

Figura 7: aplicat ,ia 2
vol <- vol+1; mass <- mass +rho(z);
volsq <- volsq+1;
masssq <- masssq +rho(z)^2;
}
}
volumeOFObject <- (vol/N)*volumeOfBox;
volumeOFObject
volvar <- (1/N)*(volsq/N-(vol/N)^2)*volumeOfBox^2;
volvar
volstd <- sqrt(volvar);
volstd
38

massOFObject <- (mass/N)*volumeOfBox;
massOFObject
massvar <- (1/N)*(masssq/N-(mass/N)^2)*volumeOfBox^2;
massvar
massstd <- sqrt(massvar);
massstd
Rezultatul afis ,at:
> volumeOFObject <- (vol/N)*volumeOfBox;
> volumeOFObject
[1] 1000
> volvar <- (1/N)*(volsq/N-(vol/N)^2)*volumeOfBox^2;
> volvar
[1] 0
> volstd <- sqrt(volvar);
> volstd
[1] 0
39

> massOFObject <- (mass/N)*volumeOfBox;
> massOFObject
[1] 2320.461
> massvar <- (1/N)*(masssq/N-(mass/N)^2)*volumeOfBox^2;
> massvar
[1] 0.8425619
> massstd <- sqrt(massvar);
> massstd
[1] 0.9179117
(3)Simularea lui e. Fie o mult ,ime de numere aleatoare iar Ns a e primul care este mai
mare dec^ at predecesorul s au. Adic a:
N=min(n:n2;Un>Un1)
Acum
PfN >ng=PfU1U2


Ung=1
n!
unde ultima egalitatea rezult a din faptul c a toate aranjamentele posibile ale lui
40

U1;


;Unsunt la fel de probabile. Cum
PfN=ng=PfN >n1gPfN >ng=1
(n1)!1
n!=n1
n!
s,i deci
M[N] =1X
n=21
(n2)!=e
Iar,
M[N2] =1X
n=2n
(n2)!=1X
n=22
(n2)!+1X
n=2n2
(n2)!= 2e+1X
n=21
(n3)!= 3e
As,adar,
D2(N) = 3ee20:7658:
N = 100000000;
n = N;
m = 0;
i = 0;
maxL = 0;
f = 0;
while n > 0
m = m + rand;
i = i + 1;
if s > 1
if i > maxL
f(i) = 1;
maxL = i;
else
f(i) = f(i) + 1;
end
41

i = 0;
s = 0;
n = n – 1;
end
end
disp ((1:maxL)*f'/sum(f))
bar(f/sum(f))
grid on
f/sum(f)
Rezultatul afisat:
>> e
2.7185
ans =
Columns 1 through 6
0 0.4999 0.3334 0.1250 0.0334 0.0069
Columns 7 through 11
0.0012 0.0002 0.0000 0.0000 0.0000
>>
42

1 2 3 4 5 6 7 8 9 10 1100.050.10.150.20.250.30.350.40.450.5Figura 8: Histograma
(4)Pentru funct ,ia:
h(x) = [cos(50x) +sin(20x)]2(5.1)
s a consider am evaluarea integralei peste [0 ;1]. Poate  v azut a ca o medie uniform a,
iar apoi vom genera U1;U2;


;Uni.i.d.U(0;1) variabile aleatoare s ,i aproxim am
R
h(x)dxcuPh(Ui)
n.
> h=function(x){(cos(50*x)+sin(20*x))^2}
> par(mar=c(2,2,2,1) ,mfrow=c(2,1))
> curve(h,xlab="Function" ,ylab="",lwd=2)
> integrate(h,0,1)
0.9652009 with absolute error < 1.9e-10
> x=h(runif(10^4))
> estint=cumsum(x)/(1:10^4)
> esterr=sqrt(cumsum((x-estint)^2))/(1:10^4)
> plot(estint,xlab="Mean and error range", type="l",lwd=+2,ylim=mean(x)+20*c(-esterr[10^4],esterr[10^4]),ylab="")
> lines(estint+2*esterr,col="gold",lwd=2)
> lines(estint-2*esterr,col="gold",lwd=2)
43

0.0 0.2 0.4 0.6 0.8 1.00 1 2 3
FunctionFigura 9: Aproximarea integralei
0 2000 4000 6000 8000 100000.8 0.9 1.0 1.1
Mean and error range
Figura 10: Aproximarea mediei cu cateva erori minime
44

(5)Algoritmi Monte-Carlo. Urm atorul exemplu ilustreaz a vericarea unor identit at ,i.
[8]Fief(x1;


;xn) un polinom cu coecient ,i rat ,ionali denvariabile de grad cel
multk^ n ecare dintre variabile. Ne dorim s a veric am dac a f0. Ideea este s a
^ nlocuim variabilele cu numere aleatoare s ,i s a calcul am valoarea polinomului. Dac a
valoarea aceasta nu este 0, polinomul nu poate  identic nul. Dac a pentru un num ar
de repat ari sucient de mare, se obt ,ine de ecare dat a valoarea 0, probabilitatea ca
polinomul s a e identic nul este mic a. Vom alege pentru variabile valori ^ ntregi din
intervalul [0 ;N1], independente s ,i uniform distribuite. Atunci avem urm atorul
rezultat:
Lema 5.1. [8](Schwartz) Dac a f nu este identic nul s ,i valorileisunt independente
s,i uniform distribuite ^ n intervalul [0;N1]atunci:
P(f(1;


;n) = 0)kn
N
Demonstrat ie: Folosim induct ,ia matematic a dup a n. Pentru n=1, lema este
adev arat a (un polinom ^ ntr-o variabil a de grad kpoate avea cel mult kr ad acini).
Fien>1 s,i ordon am f astfel:
f=f0+f1x1+f2x12+


+ftx1t
undef0;


;ftsunt polinoame^ n variabile x2;


;xn, iar termenul ftnu este identic
0 s,itk. Aplic am formula probabilit at ,ii totale s ,i avem c a:
P(f(1;


;n) = 0)P(f(1;


;n) = 0jft(1;


;n) = 0)P(ft(1;


;n) = 0)+
+P(f(1;


;n) = 0jft(1;


;n)6= 0)P(ft(1;


;n)6= 0)P(ft(1;


;n) = 0)+
+P(f(1;


;n) = 0jft(1;


;n)6= 0)
Primul termen se poate estima cu ajutorul ipoteza induct ,iei, iar al doilea este cel
45

multk
N, deoarece 1este independent a de variabilele 2;


;ns,i atunci, dac a
ultimele sunt xate astfel ^ nc^ at ft6= 0 s ,ifca polinom ^ n x1nu este identic nul.
Atunci probabilitatea ca 1s a e r ad acin a este cel multk
N. Deci
P(f(1;


;n) = 0)k(n1)
N+k
Nkn
N
Acesta ne conduce la urm atorul algoritm:
1. Calcul am f(1;


;n) pentru valori ^ ntregi 1, numere aleatoare independente, dis-
tribuite uniform discret ^ n intervalul [0 ;2K]
2. Dac a obt ,inem o valoare diferit a de 0, ne oprim: f este identic nul.
3. Dac a obt ,inem valoarea 0, repet am calculul. Dac a se obt ,ine de pre multe ori valoarea
0, ne oprim s ,i spunem c a fe identic nul.
Observat ,ie 5.2. [8] Dac a num arul de repet ari de repet arieste l, probabilitatea ca algo-
ritmul s a decid a eronat c a f0este de<2l, deoarece probabilitatea de a gres ,i la o
^ ncercare este1
2iar ^ ncerc arile sunt independente.
Vom verica identitatea
(a+b+c)3a3b3c3= 3(a+b)(a+c)(b+c)
Rutina verifident aplic a lema lui Schwarz
#verificare unor identitati – Lema lui Schwarz
verifident<-function(f,n,k) {
for (j in 1:100) {
x <- sample(0:(k*n), n, replace = TRUE);
y <- f(x);
if (y!=0) return(F);
46

}
return(T);}
Vericarea efectiv a c a
f(a;b;c ) = (a+b+c)3a3b3c33(a+b)(a+c)(b+c)0
se face cu
f2<-function(x) (x[1]+x[2]+x[3])^3-x[1]^3-x[2]^3-x[3]^3
+ -3*(x[1]+x[2])*(x[1]+x[3])*(x[2]+x[3]);
r2<-verifident(f2,3,3)
r2
[1] TRUE
Probabilitatea de eroare este mai mic a dec^ at 2100.
47

Bibliograe
[1] Christian P. Robert, George Casella, Introducing Monte Carlo Methods with R ,
SPRINGER, New York , 2010
[2] Diane O'Leary, Scientic Computing with Case Studies , SIAM, Philadelphia, 2009
[3] D.D. Stancu, Gh. Coman, P. Blaga, Analiz a numeric a s ,i Teoria Aproxim arii , Vol.
II, PRESA UNIVERSITAR A CLUJEAN A, Cluj-Napoca, 2002.
[4] James E. Gentle, Random Number Generation and Monte Carlo Methods ,
SPRINGER, New York, 2003
[5] Natalia C. Ros ,ca,Monte Carlo and Quasi-Monte Carlo Methods with Applications ,
CLUJ UNIVERSITY PRESS, Cluj-Napoca, 2009
[6] Natalia C. Ros ,ca, notit ,ele din cadrul cursului de Probabilit at ,i s,i Statistic a.
[7] N. Madras, Lectures on Monte Carlo Methods ,AMERICAN MATHEMATICAL SO-
CIETY, 2002
[8] Radu Tr^ mbit ,as,note de curs http://math.ubbcluj.ro/ ~tradu/applmstat.htm
[9] Sheldon M. Ross, Simulation, ACADEMIC PRESS , Oxford, 2012
[10] S. M. Ermakov, Metoda Monte Carlo s ,i probleme ^ nrudite , ED. TEHNIC A,
Bucures ,ti, 2009
48