Metodo Monte Carlo

Da Wikipedia, l'enciclopedia libera.
visualizzazione tridimensionale di una simulazione estremamente grande di un problema di Instabilità di Rayleigh-Taylor - Lawrence Livermore National Laboratory

Il Metodo Monte Carlo fa parte della famiglia dei metodi statistici non parametrici. È utile per superare i problemi computazionali legati ai test esatti (ad esempio i metodi basati sulla distribuzione binomiale e calcolo combinatorio, che per grandi campioni generano un numero di permutazioni eccessivo).

Il metodo è usato per trarre stime attraverso simulazioni. Si basa su un algoritmo che genera una serie di numeri tra loro incorrelati, che seguono la distribuzione di probabilità che si suppone abbia il fenomeno da indagare. L'incorrelazione tra i numeri è assicurata da un test chi quadrato.

La simulazione Monte Carlo calcola una serie di realizzazioni possibili del fenomeno in esame, con il peso proprio della probabilità di tale evenienza, cercando di esplorare in modo denso tutto lo spazio dei parametri del fenomeno. Una volta calcolato questo campione casuale, la simulazione esegue delle 'misure' delle grandezze di interesse su tale campione. La simulazione Monte Carlo è ben eseguita se il valore medio di queste misure sulle realizzazioni del sistema converge al valore vero.

Le sue origini risalgono alla metà degli anni 40 nell'ambito del Progetto Manhattan. I formalizzatori del metodo sono Enrico Fermi, John von Neumann e Stanisław Marcin Ulam[1], il nome Monte Carlo fu inventato in seguito da Nicholas Constantine Metropolis in riferimento alla nota tradizione nei giochi d'azzardo del mini stato omonimo nel sud della Francia, l'uso di tecniche basate sulla selezione di numeri casuali è citato già in un lavoro di Lord Kelvin del 1901 ed in alcuni studi di William Sealy Gosset[1].

L'algoritmo Monte Carlo è un metodo numerico che viene utilizzato per trovare le soluzioni di problemi matematici, che possono avere molte variabili e che non possono essere risolti facilmente, per esempio il calcolo integrale. L'efficienza di questo metodo aumenta rispetto agli altri metodi quando la dimensione del problema cresce.

Un primo esempio di utilizzo del metodo Monte Carlo è rappresentato dall'esperimento dell'ago di Buffon e forse il più famoso utilizzo di tale metodo è quello di Enrico Fermi, quando nel 1930 usò un metodo casuale per problemi di trasporto neutronico[1].

Descrizione generale[modifica | modifica sorgente]

Il metodo Monte Carlo può essere illustrato come una Battaglia navale. Prima un giocatore fa alcuni colpi a caso. Successivamente il giocatore applica alcuni algoritmi (es. la corazzata è di quattro punti nella direzione verticale o orizzontale). Infine, sulla base dei risultati del campionamento casuale e degli algoritmi il giocatore può determinare le posizioni probabili delle navi degli altri giocatori.

Non c'è un solo metodo Monte Carlo; il termine descrive invece una classe di approcci molto utilizzati per una larga categoria di problemi. Tuttavia, questi approcci tendono a seguire un particolare schema:

  1. Definire un dominio di possibili dati in input.
  2. Generare input casuali dal dominio con una certa distribuzione di probabilità determinate.
  3. Eseguire un calcolo deterministico utilizzando i dati in ingresso (input).
  4. Aggregare i risultati dei calcoli singoli nel risultato finale.

Integrazione[modifica | modifica sorgente]

I metodi deterministici di integrazione numerica operano considerando un numero di campioni uniformemente distribuiti. In generale, questo metodo lavora molto bene per funzioni di una variabile. Tuttavia, per funzioni di vettori, i metodi deterministici di quadratura possono essere molto inefficienti. Per integrare numericamente una funzione di un vettore bidimensionale, sono richieste griglie di punti equispaziati sulla superficie stessa. Per esempio una griglia di 10x10 richiede 100 punti. Se il vettore è a 100 dimensioni, la stessa spaziatura sulla griglia dovrebbe richiedere 10100 punti– questo potrebbe essere troppo dispendioso computazionalmente. Le 100 dimensioni non hanno un significato irragionevole, poiché in molti problemi di fisica, una "dimensione" è equivalente a un grado di libertà.

I metodi di Monte Carlo forniscono una soluzione a questo problema di crescita esponenziale del tempo. Finché la funzione in questione ha un buon comportamento, può essere valutata selezionando in modo casuale i punti in uno spazio 100-dimensionale, e prendendo alcune tipologie di medie dei valori della funzione in questi punti. Per il teorema del limite centrale, questo metodo mostrerà ordine di convergenza 1/\sqrt{N}; per esempio quadruplicando il numero dei punti equispaziati dimezza l'errore, nonostante il numero delle dimensioni.

Una caratteristica di questo metodo è quella di scegliere i punti in modo casuale, ma vengono scelti con maggior probabilità i punti che appartengono alle regioni che contribuiscono maggiormente al calcolo dell'integrale rispetto a quelli che appartengono a regioni di basso contributo. In altre parole, i punti dovrebbero essere scelti secondo una distribuzione simile in forma alla funzione integranda. Comprensibilmente, fare ciò è difficile tanto quanto risolvere l'integrale, ma ci sono altri metodi di approssimazione possibili: a partire da quelli che costruiscono una funzione integrabile simile a quella da integrare, fino ad arrivare ad una delle procedure adattive.

Un simile approccio implica l'uso di low-discrepancy sequences piuttosto del metodo quasi-Monte Carlo. I metodi Quasi-Monte Carlo spesso possono essere più efficienti come metodi di integrazione numerica poiché la successione di valori generata riempie meglio l'area e le successive valutazioni possono far convergere più velocemente la simulazione alla soluzione.


Discussione analitica[modifica | modifica sorgente]

Per un modello stocastico sia θ la quantità da determinarsi. Si esegua una simulazione, generando la variabile casuale X1 in modo che θ sia il valore atteso di X1. Consideriamo una seconda simulazione, generando una variabile casuale X2 tale che il suo valore atteso sia sempre θ. Proseguiamo con k simulazioni, generando fino a k variabili casuali Xk con E[Xk] = θ. Come stimatore di θ possiamo prendere la media aritmetica delle k variabili casuali generate, cioè

 X = \frac{\sum_{i=1}^k X_i}{k}

in quanto è ovviamente E[X] = θ. Qual è il valore più appropriato di k? Supponiamo di avere n variabili aleatorie indipendenti, X1, ..,Xn aventi la stessa distribuzione. Sia σ2 la varianza della variabile Xi e θ il valore atteso (E[Xi] = θ, Var(Xi) = σ2). La media campionaria X viene definita da

 X = \frac{\sum_{i=1}^n X_i}{n}

Il suo valore atteso è:

E[X] = \frac{\sum_{i=1}^n E[X_i]}{n} = \theta

Quindi X è uno stimatore non distorto (cioè con valore atteso uguale a quello del parametro) di θ. La sua varianza, usando la formula di Bienaymé è:

Var(X) = E[(X - \theta)^2] = Var \left [ \frac{\sum_{i=1}^n X_i}{n} \right ] \quad = \frac{\sum_{i=1}^n Var(X_i)}{n^2}= \frac{\sigma^2}{n}

Pertanto X è una variabile aleatoria con media θ e varianza σ2/n; ne segue che X è uno stimatore efficiente quando σ/√n è piccolo. Fissata una tolleranza per σ2/n ed avendo stimato σ2 si può in tal modo stimare n.

Si può imporre che il valore atteso ottenuto con lo stimatore stia dentro un ben definito intervallo di confidenza. Si può a tale scopo utilizzare una conseguenza del teorema del limite centrale. Sia X1, X2, …, Xn …, una successione di variabili casuali indipendenti e distribuite identicamente aventi la media finita μ e la varianza finita σ2. Allora

\lim_{n \to \infty}P \left ( \frac {X_1+ X_2+...+ X_n -n\mu}{\sigma \sqrt{n}} \leq x \right )  = \Phi(x)

dove Φ(x) è la funzione di distribuzione di una variabile casuale normale standard,

 \Phi(x) = \frac{1}{\sqrt{2 \pi}}\int_{-\infty}^{x}e^{ \frac{-y^2}{2}}\, dy


Quando n>>1 il teorema del limite centrale ci dice che la variabile

Z = \frac{(X - \theta)}{(\frac{\sigma}{\sqrt{n}})}

è approssimativamente distribuita come una variabile aleatoria normale unitaria, indicata con N(0,1), cioè con media zero e varianza 1. Sia ora zα, dove 0< α <1, quel numero tale che, per una variabile normale unitaria, si abbia P(Z > zα ) = α Allora, dal teorema del limite centrale si ha che , asintoticamente per n grande

P \{ X - z \left ( \frac{\alpha}{2} \right ) \frac{\sigma}{\sqrt{n}} < \theta < X + z \left ( \frac{\alpha}{2} \right ) \frac{\sigma}{\sqrt{n}} \} = 1 - \alpha

Che afferma che la probabilità che la media θ sia compresa nell'intervallo

 [X - z \left ( \frac{\alpha}{2} \right ) \frac{\sigma}{\sqrt{n}}, X + z \left ( \frac{\alpha}{2} \right ) \frac{\sigma}{\sqrt{n}} ]

è (1 - α). Perciò, assegnato 1-α e conoscendo σ, si può stimare il minimo valore di n necessario.

Nasce quindi il problema di come stimare la varianza σ2 = E[(X - θ)2]

Definizione. La varianza del campione S2 è definita da

 S^2 = \sum_{i=1}^n \frac{(X_i - X)^2}{(n - 1)}

Vale il seguente risultato.

Proposizione. E[S2]= σ2 Infatti si ha:

 \sum_{i=1}^n (X_i - X)^2 = \sum_{i=1}^n (X_i^2 - n X^2)

ne segue

 (n-1)E[S^2] = E[ \sum_{i=1}^n X_i^2 ] - n E [X^2] = n E [X_i^2] - n E [X^2]

Per una variabile aleatoria si ha:

 E[Y^2] = Var(Y) + (E[Y])^2

E quindi

E[X_i^2] = Var(X_i) + (E[Xi])^2 = \sigma^2 + \theta^2

Inoltre

 E[X^2] = Var(X) + (E[X])^2 = \frac{\sigma^2}{n} + \theta^2

Ne segue

 (n -1)E[S^2] = n \sigma^2 + n \theta^2 - \sigma^2- n \theta^2 = (n -1) \sigma^2

Supponiamo ora di avere n variabili aleatorie indipendenti X1, X2, …, Xn aventi la stessa funzione di distribuzione F e di volere stimare il parametro θ(F) (per evidenziare che tale quantità deve essere calcolata rispetto alla funzione di distribuzione F). Sia g(X1, X2, …, Xn) lo stimatore proposto per θ(F); se questo non corrisponde al valore medio, il metodo precedentemente esposto per stimare la varianza dello stimatore non si può applicare. Vediamo come si può stimare l'errore quadratico medio che si commette quando si usa questo stimatore:

 EQM(F) = E_F[(g (X_1, X_2, ... , X_n) - \theta(F))^2]

Dove il pedice F significa che il valore d'aspettazione viene calcolato rispetto alla funzione di distribuzione F che per il momento è incognita.

Un metodo per stimare tale quantità è quello del bootstrap, utilizzando la funzione di distribuzione empirica Fe(x) definita da:

 F_e(x) = \frac{(numero\; degli\; i: X_i \le x)}{n}

La legge forte dei grandi numeri afferma che per n molto grande, con probabilità 1, Fe(x) tende a F(x). Allora un valore approssimato di EQM(F) è dato da (approssimazione di bootstrap):

 EQM(F) = E_{Fe}[(g (X_1, X_2, ... , X_n) - \theta(F_e))^2]

Va rilevato, da un punto di vista operativo, che il dimensionamento della simulazione si supera facilmente grazie alla crescente disponibilità di potenza di calcolo. In altre parole, procedendo all'uso del metodo su calcolatore, sarà sufficiente generare una serie di prove di ampiezza sicuramente ridondante per assicurarsi la significatività della stima.


Esempi[modifica | modifica sorgente]

Rendimento di un titolo[modifica | modifica sorgente]

Si voglia stimare il rendimento mensile di un titolo azionario. Il titolo esiste da cinque anni, quindi si hanno a disposizione solo 60 rendimenti mensili. Supponiamo che i rendimenti si distribuiscano seguendo una variabile casuale normale.

Calcoliamo:

Con un modello di regressione lineare cercheremo di stimare la media a un mese. Successivamente, si andranno a generare attraverso l'algoritmo Monte Carlo una serie di medie "sperimentali" che saranno ricavate da una distribuzione normale (perché si è ipotizzato che i rendimenti seguano questa distribuzione) con media pari alla media stimata e scarto quadratico medio pari allo scarto quadratico medio campionario a un mese.

Una strategia per procedere e stimare la vera media del fenomeno, a questo punto, può essere quella di ricavare la media generale di tutte le medie sperimentali ottenute. I dati ottenuti forniscono stime tanto migliori quanto maggiore è il numero delle prove fatte.

Determinazione del valore π[modifica | modifica sorgente]

Monte-Carlo method pi.svg

Sia M un punto di coordinate (x,y) con 0<x<1 e 0<y<1.

Scegliamo casualmente i valori di x e y.

Sia x^2+y^2<1 allora il punto M appartiene al disco di centro (0,0) di raggio 1.

La formula per determinare l'area di un disco è il raggio elevato al quadrato per π. Nell'esempio il raggio è pari a uno e quindi l'area di interesse è 1*π = π. Il punto può cadere solo in uno dei quattro quadranti del disco e quindi la probabilità che cada all'interno del disco è π/4.

Facendo il rapporto del numero dei punti che cadono nel disco con il numero dei tiri effettuati si ottiene un'approssimazione del numero π/4 se il numero dei tiri è grande.

Eseguendo numericamente l'esempio si ottiene un andamento percentuale dell'errore mostrato nel grafico sottostante.

Andamento % dell'errore tra il pi-greco teorico e il pi-greco calcolato. Il programma ha eseguito 1370 milioni di lanci. Si noti che all'inizio l'errore è molto elevato ma rapidamente tende a decrescere. Essendo un metodo statistico ci possono essere dei temporanei innalzamenti dell'errore ma la tendenza è la sua diminuzione all'aumento dei lanci.

Determinazione della superficie di un lago[modifica | modifica sorgente]

Questo è un esempio classico della divulgazione del metodo Monte-Carlo. Sia data una zona rettangolare o quadrata di cui la lunghezza dei lati è conosciuta. Al centro di quest'area si trova un lago la cui superficie è sconosciuta. Grazie alle misure dei lati della zona, si conosce l'area del rettangolo. Per determinare l'area del lago, si chiede ad una truppa armata di tirare X colpi di cannone in modo aleatorio su questa zona. Contiamo in seguito il numero N di palle che sono restate sulla terra, possiamo quindi determinare il numero di palle che sono cadute dentro il lago: X-N. È sufficiente quindi stabilire un rapporto tra i valori:

\frac{{\it superficie}_{\,terreno}}{{\it superficie}_{\,lago}} = \frac{X}{X-N}
{\it superficie}_{\,lago} = \frac{(X-N) \cdot {\it superficie}_{\,terreno}}{X}

Per esempio, se il terreno ha superficie di 1000 m2, e supponiamo che l'armata tiri 500 palle e che 100 proiettili sono caduti dentro il lago allora la superficie del lago è di: 100*1000/500 = 200 m2.

Stima della superficie del lago grazie a dei tiri di cannone casuali

Naturalmente, la qualità della stima migliora aumentando il numero dei tiri ed assicurandosi che l'artiglieria non miri sempre lo stesso posto ma copra bene la zona. Questa ultima ipotesi coincide con l'ipotesi di avere un buon generatore di numeri aleatori, questa condizione è indispensabile per avere dei buoni risultati con il metodo Monte Carlo. Un generatore distorto è un po' come un cannone il cui tiro tende a indirizzarsi verso alcune zone rispetto ad altre: le informazioni che se ne possono ricavare sono anch'esse distorte.

Alcune applicazioni[modifica | modifica sorgente]

Il metodo è molto usato in varie discipline. Tra le possibili applicazioni: fisica statistica e ingegneria, dove si presta molto bene a risolvere problemi legati, ad esempio, alla fluidodinamica; in economia e finanza per prezzare i derivati e le opzioni non standard; in informatica, per simulare l'illuminazione naturale; in chimica computazionale il Monte Carlo quantistico è un metodo per la determinazione della struttura elettronica; ecc…

È molto potente se usato in combinazione con altri metodi non parametrici come il resampling.

Un altro esempio particolare dell'utilizzo del Metodo Monte Carlo è l'impiego del metodo nell'analisi scacchistica. Negli ultimi anni alcuni programmi scacchistici in commercio implementano delle opzioni d'analisi che utilizzano "Monte Carlo analisi". Per valutare una posizione, si fanno giocare al computer migliaia di partite partendo dalla posizione da analizzare, facendo eseguire al PC delle mosse "random" (una scelta casuale tra le mosse che il programma giudica più logiche). La media dei risultati ottenuti in queste partite è un'indicazione plausibile della mossa migliore. (Fonte: Chessbase).

Il metodo Monte Carlo trova un'applicazione di successo ancora maggiore nel gioco del go. In tal caso l'applicazione del metodo è ancora più diretta: molti programmi raggiungono un buon livello di gioco senza aver bisogno di restringere la selezione casuale a un sottoinsieme di mosse legali determinato da una funzione di valutazione.


Note[modifica | modifica sorgente]

  1. ^ a b c Carlo Jacoboni e Paolo Lugli, The Monte Carlo Method for Semiconductor Device Simulation - Springer-Verlag

Bibliografia[modifica | modifica sorgente]

  • W.K Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrikam, 1970; 57: 97-109
  • Bernd A. Berg, Markov Chain Monte Carlo Simulations and Their Statistical Analysis (With Web-Based Fortran Code), World Scientific 2004, ISBN 981-238-935-0.
  • P. Kevin MacKeown, Stochastic Simulation in Physics, 1997, ISBN 981-3083-26-3
  • Harvey Gould & Jan Tobochnik, An Introduction to Computer Simulation Methods, Part 2, Applications to Physical Systems, 1988, ISBN 0-201-16504-X
  • C.P. Robert and G. Casella. "Monte Carlo Statistical Methods" (second edition). New York: Springer-Verlag, 2004, ISBN 0-387-21239-6
  • Makers of commercial packages which implement Monte Carlo algorithms in include Palisade Corporation (@Risk), Decisioneering (Crystal Ball) and Vanguard Software (DecisionPro)
  • Mosegaard, Klaus., and Tarantola, Albert, 1995. Monte Carlo sampling of solutions to inverse problems. J. Geophys. Res., 100, B7, 12431-12447.
  • Tarantola, Albert, Inverse Problem Theory (free PDF version), Society for Industrial and Applied Mathematics, 2005. ISBN 0-89871-572-5
  • Morin, L. Richard, Monte Carlo Simulation in the Radiological Sciences, CRC Press, ISBN 0-8493-5559-1.

Voci correlate[modifica | modifica sorgente]

Altri progetti[modifica | modifica sorgente]

Collegamenti esterni[modifica | modifica sorgente]

Software Statistici[modifica | modifica sorgente]