Implementazioni di algoritmi/Metodo Monte Carlo

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 rappresentativo, 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.

Da un altro punto di vista le simulazioni Monte Carlo non sono altro che una tecnica numerica per calcolare integrali.

Le sue origini risalgono alla metà degli anni 40 all'interno del Progetto Manhattan. I formalizzatori del metodo sono John von Neumann e Stanisław Marcin Ulam, il nome Monte Carlo fu assegnato in seguito da Nicholas Constantine Metropolis in riferimento al celebre casino.

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 nel1930 usò un metodo casuale per calcolare le proprietà del neutrone.

Integrazione

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 su una superficie bi-dimensionale. Per esempio una griglia di 10x10 richiede 100 punti. Se il vettore ha dimensione 100, la stessa spaziatura sulla griglia dovrebbe richiedere 10¹⁰⁰ punti– questo potrebbe essere troppo dispendioso computazionalmente. Le 100 dimensioni non hanno significato ragionevole, poiché in molti problemi di fisica, una "dimensione" è equivalente al 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.

Esempio

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:

• Media campionaria
• Scarto quadratico medio campionario, su base giornaliera (che poi si adatterà con la formula della radice quadrata del tempo al periodo mensile.)

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.

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; in informatica, per simulare l'illuminazione naturale; ecc…

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

Discussione analitica

Per un modello stocastico sia θ la quantità da determinarsi. Si esegua una simulazione, generando la variabile casuale X₁ in modo che θ sia il valore atteso di X₁. Consideriamo una seconda simulazione, generando una variabile casuale X₂ tale che il suo valore atteso sia sempre θ. Proseguiamo con k simulazioni, generando fino a k variabili casuali X_k con E[X_k] = θ. 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, X₁, ..,X_n aventi la stessa distribuzione. Sia σ² la varianza della variabile X_i e θ il valore atteso E[X_i] = θ Var(X_i) = σ² 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 "unbiased" (cioè con valore atteso uguale a quello del parametro) di θ. La sua varianza è:

$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 σ²/n; ne segue che X è uno stimatore efficiente quando σ/√n è piccolo. Fissata una tolleranza per σ²/n ed avendo stimato σ² 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 X₁, X₂, …, X_n …, una successione di variabili casuali indipendenti e distribuite identicamente aventi la media finita μ e la varianza finita σ². Allora

$\lim _{n\to \infty }P\left({\frac {X_{1}+X_{2}+...+X_{n}-n\mu }{\sigma {\sqrt {n}}}}\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 σ² = E[(X - θ)²]

Definizione. La varianza del campione S² e’ definita da

$S^{2}=\sum _{i=1}^{n}{\frac {(X_{i}-X)^{2}}{(n-1)}}$

Vale il seguente risultato.

Proposizione. E[S^{2]= σ2 Infatti si ha:}

$\sum _{i=1}^{n}(X_{i}-X)^{2}=\sum _{i=1}^{n}(X_{i}^{2}-nX^{2})$

ne segue

$(n-1)E[S^{2}]=E[\sum _{i=1}^{n}X_{i}^{2}]-nE[X^{2}]=nE[X_{i}^{2}]-nE[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 X₁, X₂, …, X_n 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(X₁, X₂, …, X_n) 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 F_e(x) definita da:

$F_{e}(x)={\frac {(numero\;degli\;i:X_{i}\leq x)}{n}}$

La legge forte dei grandi numeri afferma che per n molto grande, con probabilità 1, F_e(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.

Esempio: determinare il valore π

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

Scegliamo aleatoriamente 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 probabilità che il punto M appartenga al 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.

Esempio: determinare la superficie di un lago

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 {superficie_{~terreno}}{superficie_{~lago}}}={\frac {X}{X-N}}$

$superficie_{~lago}={\frac {(X-N)\cdot superficie_{~terreno}}{X}}$

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

Estimation de la surface du lac grâce à des tirs d'artillerie aléatoires

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 è come un cannone che tira sempre nello stesso punto: le informazioni che genera sono ridotte.

Bibliografia

• 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 0898715725

Collegamenti Internet

• http://epicws.epm.ornl.gov/rsic.html ORNL Radiation Safety Information Computational Center (RSICC)
• http://www.nea.fr/html/dbprog/ NEA Data Bank Computer Program Services
• http://homepages.nyu.edu/~sl1544/articles.html Simon Leger, Article on Monte Carlo techniques applied to finance
• http://homepages.nyu.edu/~sl1544/MonteCarloNuls.pdf Introduction aux techniques de Monte Carlo appliquees a la finance et introduction aux techniques plus avancees, Simon Leger, website at http://homepages.nyu.edu/~sl1544/

Software Statistici

• R packages[1]:MCMCpack Package