Metoda Monte Carlo

Metoda Monte Carlo

Introducere in

simularea Monte Carlo

O tehnica ce a avut un impact mare in diferite domenii ale stiintei computationale este tehnica numita "simularea Monte Carlo". Aceasta tehnica isi deriva numele de la cazinourile din Monte Carlo - o simulare Monte Carlo foloseste numere aleatoare pentru a modela un fel de proces. Aceasta tehnica functioneaza bine in special cand procesul este unul in care probabilitatile care stau la baza acestuia sunt cunoscute dar in care rezultatul este mai greu de determinat. O mare parte din timpul procesoarelor ale unora dintre cele mai rapide computere din lume este folosit pentru a executa simulari Monte Carlo deoarece putem scrie unele din cele mai fundamentale legi ale fizicii dar nu putem rezolva analitic pentru probleme de interes.

Un exemplu al modului de functionare a simularii Monte Carlo in viata de zi cu zi este un proiect pe care un student l-a facut la un curs de FORTRAN - acesta a vrut sa afle cea mai buna strategie de a castiga bani la blackjack. Abordarea conventionala (folosind doar statistici) ar fi aceea de a scrie probabilitatea de a avea o anumita combinatie (exemplu ar fi sa ai un as si un cinciar si crupierul sa arate un Joker) si apoi calcularea platii asteptate pentru fiecare scenariu posibil (inca sau nicio carte in plus, dar apoi va trebui sa calculezi ce sa faci daca primesti un septar). Daca stai sa te gandesti la toate variantele, in curand va deveni coplesitor. Aceasta problema, totusi, functioneaza foarte bine ca o simulare Monte Carlo. Stim probabilitatile de baza (sa apara o anumita carte are probabilitatea de 1/52 sau daca sunt 5 pachete de carti aceasta este de 1/5*52) asa ca tot ce ne trebuie sunt "regulile" pe care sa le folosim. In acest caz, studentul a scris un program care va genera aleator un joc cu 5 pachete de carti. Apoi va "imparti" cartile de joc intre el si "crupierul". Acesta din urma respecta regulile de baza (mai cere o carte la 16 si sta la 17) si studentul a programat strategia de pariu pe care vroia sa o testeze (pe ce fel de combinatii va mai cere o carte sau va sta, va dubla, va imparti, etc.). Apoi va rula programul si-l va pune sa genereze cateva sute de perechi de carti si sa urmareasca de cate ori castiga (sau pierde) si sa afiseze rezultatul la sfarsit (asta ar dura aproape o ora pe un PC). Se pot testa strategii variate si sa se observe cum vor actiona pe termen lung.

Sa luam un exemplu simplu a unei simulari Monte Carlo pentru a ilustra tehnica. Mai intai, sa luam in considerare urmatoarea problema, vrem sa facem o simulare care ne va permite sa aflam valoarea lui Pi. Vom face acest lucru astfel: se da un patrat care are un colt in originea sistemului de coordonate si laturile au lungimea 1 - evident va avea aria 1. Acum sa luam un sfert de cerc cu raza de 1 care sa fie inscris in patrat - stim ca aria este Pi/4. Putem folosi o simulare Monte Carlo pentru a afla aria relativa a cercului si a patratului si sa inmultim aria cercului cu 4 pentru a-l afla pe Pi. In particular, modul in care vom afla aria cercului este acesta: pentru ca un punct de coordonate (X, Y) sa fie inauntrul unui cerc de raza 1, distanta acestuia fata de origine (X² + Y²) va fi mai mica sau egala cu 1. Putem genera mii de perechi aleatoare de coordonate (X, Y) si sa determinam daca fiecare dintre ele sunt inauntrul cercului. De fiecare data cand se afla inauntrul cercului, vom adauga 1 contorului. Dupa ce am generat un numar mare de puncte, ratia numerelor de puncte care se afla in cerc / numarul total de puncte generate se va apropia de ratia dintre aria cercului / aria patratului. Astfel valoarea lui Pi va deveni simplu

Pi este aproximativ 4 * (Numarul de puncte din interiorul cercului) / (Numarul total de puncte generate)

Astfel putem gasi o aproximare la pentru Pi folosind matematica elementara.

Programul urmator poate fi folosit pentru a gasi o aproximare pentru Pi

% Program Matlab de cautare a lui Pi folosind numere aleatoare

% Tom Huber, 15 Iunie, 1996

Nrand = input('Cate numere aleatoare ');

NInside = 0;

for nloops=1:Nrand

Xrand = rand; % Genereaza un punct XY aleator

Yrand = rand;

Rrand = Xrand^2 + Yrand^2; % Gaseste originea fata de distanta

if (Rrand <= 1)

NInside = NInside + 1;

end

end

disp(['Total Generate: ' num2str(Nrand) ' Puncte interioare: '

num2str(NInside)]);

piapprox = 4*NInside/Nrand;

disp([' Pi aproximativ = ' num2str(piapprox)]);

Sa se ruleze programul cu aproape 1000 de puncte aleatorii. Cat de buna este aproximarea? Ne da acelasi rezultat de fiecare data cand rulezi programul? Rezultatul se imbunatateste daca rulam programul cu 5000 de puncte? Cat timp in plus ii ia sa calculeze?

Ne putem imbunatati semnificativ viteza programului optimizandu-l pentru Matlab. In particular, Matlab este foarte rapid in lucrul cu vectori si cu matrici de numere. In editia pentru studenti a Matlab-ului, cel mai mare vector care poate fi folosit poate avea 8192 de elemente, deci vom pune programul sa genereze 8192 de numere aleatoare cu un singur apel al functiei rand. Apoi vom determina raza pentru toate acestea simultan, folosind comanda

Rrand = Xrand.^2 + Yrand.^2;

A se observa utilizarea operatorului - acesta ridica la patrat fiecare element al vectorului separat in loc sa faca o inmultire in matrice. in final, putem verifica toate elementele vectorului pentru a vedea daca raza este < 1 folosind o singura comanda

CheckValue = Rrand<=1.;

Aceasta va crea un vector CheckValue care va avea cate un 1 pentru fiecare element care indeplineste conditia (Rrand<=1.) si 0 pentru fiecare element care nu indeplineste conditia. La sfarsit, putem determina numarul din interior prin adunarea tutuor valorilor din CheckValue.

NInside = NInside + sum(CheckValue);

Urmatorul program poate fi folosit pentru a aproxima valoarea lui Pi si va fi mult mai rapid decat versiunea precedenta.

% Program Matlab optimizat pentru a-l gasi Pi folosind numere aleatoare

% Tom Huber, 15 Iunie, 1996

Nrand = 8192; % Cea mai mare marime a unui vector in versiunea Student

Nmax = input('Cate bucle (din 8192 de numere fiecare)');

NTrand = 0;

NInside = 0;

for nloops=1:Nmax

Xrand = rand(1,Nrand); % Genereaza 8192 de puncte XY aleatoare

Yrand = rand(1,Nrand);

Rrand = Xrand.^2 + Yrand.^2; % Gaseste raza pentru fiecare din cele %8192 de puncte

CheckValue = Rrand<=1.; % Are 1 daca True & 0 daca False pentru %fiecare element

NInside = NInside + sum(CheckValue); % Totalul numerelor interioare

NTrand = NTrand + Nrand; % Totalul de perechi de numere %generate

end

disp(['Total Generate: ' num2str(NTrand) ' Pct interioare: '

num2str(NInside)]);

piapprox = 4*NInside/NTrand;

pierror = 4*sqrt(NInside)/NTrand;

disp([' Approximation to pi = ' num2str(piapprox)

' Cu eroarea ' num2str(pierror)]);

Rulati o data programul de deasupra si observati cat de rapid este decat sa rulam programul anterior de 8192 de ori. Acum rulati programul de 10 sau de 100 de ori si observati cum acuratetea aproximarii se imbunatateste.

Pentru simularile Monte Carlo, procesele sunt aleatoare, asa ca de fiecare data cand este rulat va veni cu rezultate usor diferite. Poate fi aratat ca eroarea dintr-un numar aleator de numarari generate de o simulare Monte Carlo este aproximativ radicalul numarului insusi.

Astfel, in acest caz incertitudinea din valoarea noastra a lui Pi este

pierror = 4*sqrt(NInside)/NTrand;

Aceasta arata o problema la o metoda Monte Carlo - pentru a imbunatati acuratetea rezultatelor cu un factor de 10 trebuie sa rulam intr-un timp de aproximativ 100 de ori mai mare. Concluzia este ca desi aceasta metoda de a-l gasi pe Pi este conceptual foarte usoara, nu ar fi o metoda eficienta de a-l gasi pe Pi pentru un numar mare de cifre (sunt programe care au calculat valoarea lui Pi pentru peste un milion de cifre, dar nu folosind o abordare de tip Monte Carlo deoarece i-ar fi trebuit mult prea mult timp de procesare).

Acum, pentru problema pe care o studiem, si anume determinarea climei globale, sunt cateva locuri unde o simulare Monte Carlo poate fi de folos. In particular, ne va trebui pentru a ne ajuta sa determinam temperatura medie globala si cantitatea de lumina solara care cade pe fiecare latitudine. Pentru a determina temperatura medie globala, vrem media temperaturii fiecarei latitudini, dar este evident mult mai mult pamant in regiunea de la ecuator catre o latitudine de 10° decat in regiunea de la 80° catre Polul Nord. Asadar pentru a determina temperatura medie, vrem sa analizam temperatura fiecarei benzi raportate la suprafata pamantului din acea banda. Putem face asta analitic prin integrare, dar aceasta poate fi facuta deasemenea cu o metoda Monte Carlo.

Vom modifica programul de dinainte pentru a genera puncte aleatoare (XYZ) intr-un cub cu fetele de o unitate. Apoi vom determina daca acestea sunt pe suprafata unei sfere de raza 1 folosind urmatoarele:

Rrand = Xrand.^2 + Yrand.^2 + Zrand.^2;

CheckValue = Rrand<=1.01 & Rrand>=.99;

aceasta va determina daca punctele sunt pe suprafata sferei. Apoi, vom verifica daca sunt puncte in interiorul fiecarei latitudini dar si daca sunt pe suprafata sferei. Programul va incrementa un contor pentru fiecare punct care indeplinesc aceste criterii. La sfarsit, putem imparti numarul din fiecare latitudine cu numarul total de puncte care se gaseau pe suprafata pentru a gasi aria din fiecare banda.

Urmatorul program va face asta

% Program ce determina raportul de arie in benzi de latitudine pe o sfera

% Tom Huber, 25 Iunie, 1996

Theta1 = 0;

Theta2 = 90;

NSubDiv = 9; % Noua subdiviziuni a cate 10 grade fiecare

dTh = (Theta2-Theta1)/NSubDiv; % Latimea fiecarei diviziuni (10 Grade)

ThLow = Theta1:dTh:Theta2-dTh; % Limita inferioara pt fiecare regiune %(0,10,20..80)

ThHigh = Theta1+dTh:dTh:Theta2; % Limita superioara pt fiecare regiune %(10,20,30..90)

Nrand = 8192; % Numarul de puncte pt editia Student

Nmax = input('Cate cicluri de 8192 de valori fiecare ');

NTrand = 0;  % Initializeaza numarul total de puncte generate

NGoodPts = 0;  % Initializeaza numarul total de puncte in sfera

NZone = zeros(1,NSubDiv); % Initializeaza numarul in fiecare zona

T0 = clock;  % Urmareste timpul de procesare (pt referinta)

for nloops=1:Nmax

Xrand = rand(1,Nrand); % Genereaza puncte XYZ in spatiu

Yrand = rand(1,Nrand);

Zrand = rand(1,Nrand);

Rrand = Xrand.^2 + Yrand.^2 + Zrand.^2; % Distanta fata de origine

CheckValue = Rrand<=1.01 & Rrand>=.99; % Este pe suprafata sferei

NGoodPts = NGoodPts + sum(CheckValue); % Totalul pe suprafata

Lat = asin(Zrand)*180/pi; % Latitudinea fiecarui pct

for i=1:NSubDiv % Treci prin toate pctele

NZoneCheck = Lat < ThHigh(i) & Lat >= ThLow(i); % Verifica daca in %latitudine

NZoneCheck = NZoneCheck .* CheckValue; % sau pe suprafata

NZone(i) = NZone(i) + sum(NZoneCheck); % Da, adauga la suma

end

NTrand = NTrand + Nrand; % Numarul total de puncte %generate

end

T0 = clock - T0;  % Timpul procesorului

disp(['Total Generate: ' num2str(NTrand) ' Pct bune: '

num2str(NGoodPts) ' Secunde: ' num2str(T0)]);

fLatitude = NZone/NGoodPts;

fError = fLatitude./sqrt(NZone);

fActual = sin(ThHigh*pi/180.)-sin(ThLow*pi/180.);

disp('Sumar pentru zone');

disp('Unghi de jos, Unghi de sus, Raport din banda simulat, Incertitudine,');

disp(' Functie actuala (folosind calculul)');

disp([ ThLow' ThHigh' fLatitude' fError' fActual']);