momente şi schiţe de informatică şi matematică
To attain knowledge, write. To attain wisdom, rewrite.

Aplicaţii elementare ale metodei Monte Carlo

C++ | Python
2013 dec

[1] A. M. Iaglom, I. M. Iaglom - "Probleme neelementare tratate elementar", Ed. Tehnică 1983

O primă ilustrare (clasică) pentru metoda Monte Carlo constă în estimarea experimentală a numărului $\pi$ - se produc aleatoriu puncte $(x,\,y) \in [0,\,1]\times [0,\,1]$ şi se contorizează câte dintre acestea "intră" în discul $x^2+y^2 \leq 1$, ajungând astfel la o valoare aproximativă a raportului dintre aria sfertului de disc şi aria pătratului de latură 1 care îl circumscrie:

Următorul program Python constituie o listă de puncte aleatoare, din care separă apoi pe cele interioare cercului - producând rezultate precum cel redat mai sus:

import matplotlib.pyplot as plt
from pylab import *
import random as rnd
rnd.seed() # iniţializează generatorul de numere pseudo-aleatoare
samples = 5000 # numărul de puncte "aruncate" în pătratul $[0,1]\times[0,1]$
points = [(rnd.random(), rnd.random()) # listă de puncte aleatoare $(x,y)\in[0,1]\times[0,1]$
          for i in range(samples)] 
points_in = [p for p in points # separă punctele interioare cercului $x^2+y^2=1$
               if sqrt(p[0]*p[0] + p[1]*p[1]) <= 1.0]
points_ext = [p for p in points if p not in points_in]
# plotează punctele interioare cercului (magenta) şi pe cele exterioare
plt.plot([p[0] for p in points_in], [p[1] for p in points_in], 'mo',
         [p[0] for p in points_ext], [p[1] for p in points_ext], 'yo')
p_in = len(points_in) # raportează numărul de puncte "magenta" la totalul punctelor $(\approx\pi/4)$
title(str(p_in) + ' din ' + str(samples) + " puncte = " + str(4.0*p_in/samples))
angle = linspace(0, pi/2, 100) # trasează (cu cyan) sfertul de cerc
plot(cos(angle), sin(angle), 'c', linewidth=3)
plt.show() # rezultă o figură precum cea reprodusă mai sus

Este uşor (dacă ignorăm imaginea grafică) să modelăm experimentul de mai sus în C++:

#include <iostream>
#include <cstdlib>  // srand(), rand(), MAX_RAND
#include <ctime>    // time()
#include <cmath>    // sqrt()
using namespace std;

double mc_pi(int samples) {
    double success = 0; int n = samples;
    while(n --) {
        double x = static_cast<double>(rand()) / RAND_MAX;
        double y = static_cast<double>(rand()) / RAND_MAX;
        if(sqrt(x*x + y*y) <= 1)
            success += 1; // contorizează punctele interioare discului
    }
    return (4*success / samples); // aproximare a numărului $\pi$
}

int main() {
    srand(time(0));
    for(int i=1; i < 6; i++)
        cout << mc_pi(i*1000000) << '\n';
}
vb@vb:~/clasa11$ ./test
3.14039
3.14249
3.14137
3.14092
3.14227

În această manieră pot fi tratate (admiţând aproximări) şi probleme dificile de teoria probabilităţilor. De exemplu (problema 95, din [1]): considerând triunghiurile înscrise într-un acelaşi cerc, care este probabilitatea ca alegând unul la întâmplare, acesta să fie triunghi ascuţitunghic?. Ar fi de formulat un program care să simuleze de 10000 de ori (sau mai mult) alegerea aleatorie a trei valori distincte $\theta\in[0,\,2\pi)$ - pentru trei puncte $(\cos\theta,\,\sin\theta)$ ale cercului - şi să contorizeze cazurile când triunghiul format este ascuţitunghic.

Confidenţa rezultatelor (şi succesul metodei Monte Carlo) depinde esenţial de calitatea generatorului de numere aleatoare implicat; în principiu trebuie asigurată o distribuţie uniformă, având o aceeaşi probabilitate de alegere pentru fiecare valoare din gama respectivă.

Un exemplu din domeniul inteligenţei artificiale

Funcţia redată mai jos determină un răspuns "onorabil" în cursul unui joc de Hex (dar poate fi vorba şi de alte jocuri), folosind deasemenea metoda Monte Carlo: se constituie lista mutărilor posibile în poziţia curentă şi apoi, pentru fiecare mutare M dintre acestea se simulează de un anumit număr de ori continuarea jocului până la finalizarea acestuia - evaluând mutarea M prin numărul de câştiguri proprii obţinute în simulările respective; mutarea de răspuns va fi acea mutare M dintre cele astfel încercate, care obţine scorul maxim.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
const int MC_SIM = 1000; // simulări ale jocului, pentru fiecare mutare posibilă

void HexGame::my_move() { // calculatorul este la mutare
    vector<int> moves;
    gen_moves(moves); // generează răspunsurile posibile (pentru partea "Computer")
    int best = 0; // cel mai bun scor, simulând jocul după fiecare mutare posibilă
    int rsp = moves[0]; // un răspuns (indiferent), dacă în final 'best'==0 ("pierde" jocul)  
    int rest = moves.size() - 1; // câmpuri de completat alternativ, după prima mutare
    Occupy my; // harta câmpurilor ocupate de "Computer" în parcursul următor
    for(int crt=0; crt <= rest; ++crt) {
        int move = moves[crt]; // încearcă fiecare dintre mutările posibile
        make_move(move);
        if(get_end()) { // încheie ("Computer win"), dacă mutarea curentă deja câştigă jocul
            undo_move(move);
            make_move(move /_dim, move %_dim);
            return;
        }
        vector<int> cand;
        gen_moves(cand); // răspusurile posibile ale oponentului
        int n_wins = 0; // de câte ori câştigă "Computer", în cele MC_SIM simulări
        for(int t=0; t < MC_SIM; ++t) {
            shuffle(cand.begin(), cand.end(), MT); // <random> C++11; <algorithm>
            my = _player[_tomove]; // HexGame::_tomove rămâne setat pe "Computer"!
            for(int i=0; i < rest; i+=2) // "Computer" face jumătate dintre mutările rămase
                my.set(cand[i]);
            if(my_win(my)) // contorizează poziţiile finale în care "Computer" câştigă
                n_wins++;
        }
        if(n_wins > best) { // Dacă simulările pentru mutarea curentă îmbunătăţesc scorul 
            best = n_wins;  // din precedentele simulări - actualizează 'best' şi 'rsp'
            rsp = move;
        }
        undo_move(move);    // Revine la poziţia iniţială, reluând apoi simulările
    }                       // pentru următoarea mutare posibilă
    if(best==0) cout << "I know that You win!\n"; // toate simulările au avut scorul 0
    make_move(rsp /_dim, rsp %_dim); // execută mutarea cu scorul 'best'
}

_player[] este un tablou intern (membru privat al clasei HexGame) care păstrează poziţiile ocupate de fiecare jucător, pe parcursul desfăşurării jocului; make_move(move) din linia 12 înscrie mutarea 'move' în câmpul _player[_tomove], dar fără a comuta valoarea _tomove (instituită intern pentru a indica jucătorul care urmează la mutare) - în timp ce make_move(linie, coloană) din linia 36 face efectiv mutarea, comutând _tomove şi actualizând diagrama curentă a jocului.

În linia 19 se constituie lista câmpurilor rămase libere, după înscrierea mutării curente 'move'; continuarea jocului (până la finalizarea lui) echivalează - pentru cazul jocului Hex - cu ocuparea de către partea indicată de _tomove a jumătate dintre câmpurile libere (în mod implicit, celelalte sunt ocupate de către oponent). Fiecare simulare a continuării jocului constă în amestecarea aleatorie a câmpurilor libere - în linia 22, folosind funcţia shuffle() din <algorithm>, sub C++11 - şi apoi, în liniile 23-25 se completează (în variabila locală "my") poziţia curentă _player[_tomove] cu jumătate dintre câmpurile libere (parcurgând din doi în doi, lista amestecată aleatoriu a acestora).

În principiu, programul subiacent HEX Monte Carlo player este cât se poate de modest, ideea fiind că poate "să joace" acceptabil (nu-i uşor - ca jucător neexperimentat - să-l baţi) în pofida faptului că nu instrumentează de loc elemente de strategie specifice jocului (şi nu implică un mecanism minimax) - programul ştiind doar ce este o mutare legală şi cum să recunoască dacă o parte a câştigat, mutările sale fiind ghidate (fără artificii) de scorurile simulărilor Monte Carlo efectuate.

vezi Cărţile mele (de programare)

docerpro | Prev | Next