Lezione 10 e 11 : metodi Montecarlo

In questa lezione proveremo a familiarizzare con l'utilizzo di tecniche numeriche basate su numeri casuali. Per prima cosa impareremo a costruire un generatore di numeri casuali. Lo utilizzeremo quindi per generare numeri che seguano diverse distribuzioni di probabilità (uniforme, esponenziale e Gaussiana). Come applicazione utilizzeremo la generazione di numeri casuali per calcolare numericamente integrali mono e multi-dimensionali. Nella prossima lezione poi utilizzeremo un generatore di numeri casuali per simulare il comportamento di un apparato sperimentale.

ESERCIZIO 10.0 - Generatore di numeri casuali (da consegnare):

In questo esercizio proveremo a costruire un generatore di numeri casuali e a studiarne il comportamento. Scriveremo un programma che produca quattro istogrammi contenenti ciascuno 10000 numeri pseudo-casuali estratti da :

una distribuzione uniforme tra 5 e 10.
una distribuzione esponenziale tra 0 e $+\infty$ con costante $\lambda$ =1.
una distribuzione gaussiana centrata in 1 e larghezza 1 con il metodo di Box-Muller.
una distribuzione gaussiana centrata in 1 e larghezza 1 con il metodo accept-reject.

Per risolvere questo esercizio si può seguire lo schema seguente:

Scrivere una classe opportuna RandomGen per la generazione di numeri casuali. La classe dovrà avere un costruttore che accetti un seed di input e si faccia carico di inizializzare i parametri del generatore ai valori nominali. La classe dovrà inoltre contenere un metodo che implementi un generatore lineare congruenziale di base e tutti i metodi necessari per le distribuzioni richieste.
Preparare un semplice main per provare le funzionalità della classe RandomGen producendo quattro istogrammi per le distribuzioni indicate sopra.

Header file della classe RandomGen::

#ifndef __RANDOMGEN_H__
#define __RANDOMGEN_H__

class RandomGen {

public:

  RandomGen(unsigned int);

  void SetA(unsigned int a) {m_a=a;}
  void SetC(unsigned int c) {m_c=c;}
  void SetM(unsigned int m) {m_m=m;}

  double Rand( );                                  // distribuzione uniforme tra 0 e 1
  double Unif(double xmin, double xmax);           // distribuzione uniforme tra xmin e xmax
  double Exp(double lambda);                       // distribuzione esponenziale con costante lambda
  double Gaus(double mean, double sigma);          // distribuzione gaussiana (Box-Muller)
  double GausAR(double mean, double sigma, ... );  // distribuzione gaussiana (Accept-Reject)

 private:

  unsigned int m_a, m_c, m_m;
  unsigned int m_seed;

};

#endif

Il costruttore deve accettare un unsigned int come seed di input e farsi carico di inizializzare i parametri del generatore ai valori nominali

m_a = 1664525
m_c = 1013904223
m_m = $2^{31}$ .

Main per il test del generatore RandomGen:

Un possibile main per provare il generatore appena costruito è riportato di seguito :

#include "TApplication.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TAxis.h"

#include <iostream>

#include "RandomGen.h" 

int main() {

  TApplication app("app",0,0);

  RandomGen myGen(1);

  int nmax = 10000;

  TH1F unif("Uniforme","Uniforme",70,4,11) ;

  for ( int k = 0 ; k < nmax ; k++ ) {

    unif.Fill( myGen.Unif(5,10)  ) ;  

    // aggiungere le altre distribuzioni

  }

  TCanvas can2("Uniforme","Uniforme") ;
  can2.cd();
  unif.GetXaxis()->SetTitle("x [AU]");
  unif.GetYaxis()->SetTitle("N");
  unif.Draw();

  app.Run();

}

Alcuni suggerimenti qui sotto :

Generatore Lineare congruenziale

I generatori lineari congruenziali generano una sequenza di interi pseudo-casuali utilizzando la sequente formula: $n_{i} = mod ( a \, n_{i-1} + c, m)$ Quindi è necessario tenere in memoria il valore del numero intero generato al passaggio precedente. Alla prima iterazione $n_{0} = seed$

Per ottenere un numero floating point uniformemente distribuito tra 0 e 1 (con 1 escluso) è sufficiente richiedere: $d_{i} = n_{i}/m$

Ricordiamo che in C++ l'operazione mod può essere eseguita con l'operatore % che ritorna il resto della divisione (intera) per cui, per esempio, $7 \% 5 = 2$ .

Generatore esponenziale

Per generare numeri pseudo-casuali secondo la densità di probabilità esponenziale $p(x) = \begin{cases} \lambda e^{- \lambda x} & (x\ge 0) \\ 0 & (x<0) \end{cases}$ si può facilmente invertire la funzione cumulativa di $p(x)$ (metodo della funzione inversa). A partire da un numero $u$ estratto secondo la distribuzione uniforme tra 0 e 1, il numero $x = - \frac{1}{\lambda} ln (1-u)$ è distribuito proprio come la distribuzione esponenziale $p(x)$ ;

Metodo di Box Muller

Il metodo di Box-Muller ci permette di estrarre numeri secondo una densità di probabilità gaussiana sfruttando il metodo della funzione inversa bidimensionale. Siano $x$ e $y$ due variabili indipendenti distribuite normalmente (Gaussiana con media 0 e sigma 1). L'espressione della loro distribuzione di probabilità (PDF) in due dimensioni è $P(x,y) = \frac{1}{2} e^ { \frac{-(x^{2}+y^{2})}{2}}$ passando alle coordinate polari $x = r cos(\theta)$ $y = r sin(\theta)$ si ha $P(r,\theta) = \frac{1}{2} e^{ \frac{ -r^{2}}{2}}$ se calcoliamo l'integrale di tale PDF per $r \in [0,R]$ e $\theta \in [0,\Theta]$ otteniamo $F(R,\Theta) = (1 - e^{ \frac{- R^{2}} {2} }) \frac{ \Theta } {\pi}$ Definendo s e t come due variabili casuali distribuite uniformemente in [0,1], abbiamo $R = \sqrt{(-2ln(1-s)}$ $\Theta=2 \pi t$ per cui due variabili x e y distribuite normalmente possono essere generate a partire da una coppia s e t distribuita uniformemente in [0,1], secondo la formula $x = \sqrt{-2ln(1-s) } * cos(2 \pi t)$ $y = \sqrt{-2ln(1-s) } * sin(2 \pi t)$ Se ci limitiamo per comodità a considerare solo una delle due possiamo generalizzare al caso di variabile $x$ distribuita gaussianamente con media $\mu$ e larghezza $\sigma$ $x = \mu + \sigma * \sqrt{-2ln(1-s) } * cos(2 \pi t)$ .

Una possibile implementazione di tale metodo:

double RandomGen::Gaus(double mean, double sigma) {
  double s=Rand();
  double t=Rand();
  double x=sqrt(-2.*log(1.-s))*cos(2.*M_PI*t);
  return mean+x*sigma; 
}

Si noti che tecnicamente dati due numeri s e t possiamo generare due numeri x e y distribuiti gaussianamente. Nel nostro esempio ne consideriamo uno soltanto per comodità di scrittura codice, non certo per efficienza !

Metodo Accept-Reject

Il metodo accept-reject può essere utilizzato per generare numeri casuali distribuiti secondo una qualsivoglia forma funzionale. Consideriamo di voler generare numeri nell'intervallo $[a,b]$ secondo la distribuzione $f(x)$ in figura

Il metodo si basa sulla generazione di una coppia di numeri $x \in [a,b]$ e $y \in [0,max]$ dove $max$ è un numero maggiore del massimo valore assunto da $f(x)$ nell'intervallo $[a,b]$ . La coppia $x,y$ può essere facilmente generata a partire da due numeri $s$ e $t$ generati uniformemente in $[0,1]$ usando le formule $x = a + (b-a) \cdot s$ $y = max \cdot t$ Generata la coppia $x,y$ si valuta quindi $f(x)$ e si accetta $x$ se $y<f(x)$ , altrimenti si ripete la procedura. Così facendo si avrà un maggior numero di punti generati laddove $f(x)$ assume valori più grandi. Potete pensare di scrivere un metodo accept-reject specifico per la gaussiana o pensarene uno più generale che per esempio accetti una generica FunzioneBase come input.

ESERCIZIO 10.1 - Verifica del Teorema del Limite Centrale (da consegnare):

Una applicazione molto interessante dell'estrazione di numeri casuali è la verifica del teorema del limite centrale. Per fare questo possiamo generare una serie di numeri casuali uniformemente distribuiti in [0,1] e calcolare la somma eseguita su un numero N di elementi consecutivi della serie generata. Calcolare la varianza della serie di numeri generata e della serie delle somme. Verificare come questa scala con N. Si possono pensare a due versioni di questo esercizio:

semplice : passare da riga di comando il numero N di elementi da sommare. Creare due istogrammi che contengano la distribuzione dei numeri generati e la distribuzione delle somme di N elementi. Verificare come cambia la distribuzione delle somme al variare di N provando ad eseguire il codice varie volte con diversi valori di N.
completo : il codice prepara le distribuzioni (istogrammi) delle somme di N elementi con N che va da 1 a 12. Per ogni N calcola la varianza della distribuzione e genera un plot finale della varianza in funzione di N.

In entrambi i casi un numero di entries maggiore di 10000 dovrebbe andare bene.

Il Teorema del Limite Centrale

I teoremi del limite centrale sono una famiglia di teoremi di convergenza debole nell'ambito delle teoria delle probabilità. Per tutti vale l'affermazione che la distribuzione di probabilità della somma (normalizzata) di un gran numero di variabili casuali tende ad una data distribuzione regolare (attrattore), che di norma è la Gaussiana o la Lorenziana.

Nel nostro caso, verificheremo che la somma di N variabili aleatorie $x_{i}$ indipendenti e identicamente distribuite con media $\mu$ e varianza $\sigma^{2}< + \infty$ , tende a distribuirsi come una variabile casuale gaussiana con media pari a N $\mu$ e varianza N $\sigma^{2}$ , al tendere di N a infinito.

Risultati attesi

L'output del programma potrebbe essere il seguente : per prima cosa un grafico che mostri le distribuzioni delle somme e un grafico che mostri la varianza delle distribuzioni precedenti in funzione del numero di elementi sommati

ESERCIZIO 10.2 - Calcolo di integrali con metodi MonteCarlo (da consegnare):

Proviamo in questo esercizio a studiare il comportamento delle tecniche MonteCarlo per il calcolo numerico di un integrale mono-dimensionale.

Calcolare 10000 volte il valore dell'integrale di sin(x) tra [0,π] utilizzando il metodo della media con N=100 punti e fare un grafico ( istogramma ) della distribuzione dei valori ottenuti.
Estendere il punto precedente calcolando 10000 volte il valore dell'integrale di sin(x) tra [0,π] utilizzando il metodo della media a N punti con N pari a 500, 1000, 5000, 10000, 50000 e 100000 punti. Per ogni valore di N produrre il grafico della distribuzione dei 10000 valori ottenuti.
Stimare l'errore sul calcolo dell'integrale a 500, 1000, 5000, 10000, 50000, 100000 punti come deviazione standard dei 10000 valori calcolati per ogni N. Far un grafico di questo errore in funzione di N.
Assumendo che l'andamento dell'errore sia noto ( del tipo $k/ \sqrt{N}$ ) si determini quanti punti sono necessari per ottenere una precisione di 0.001. Si ripeta lo stesso lavoro con il metodo hit-or-miss.

Alcune osservazioni per lo svolgimento dell'esercizio :

Per il calcolo di integrali con metodi MonteCarlo si può decidere di scrivere una nuova classe dedicata o estendere la classe Integrale della lezione 7. Notate che in questo caso la nuova classe avrà avere un generatore di numeri casuali come data membro.
poichè il calcolo degli integrali con N molto elevato potrebbe richiedere un certo tempo, diventa molto utile spezzare il main in due main separati:
il primo calcola gli integrali e li salva in diversi files di testo (uno per ogni valore di N).
il secondo legge i files scritti dal primo e li analizza.

In questo modo possiamo agevolemnte ripetere molto volte l'analisi dei dati senza dover ricalcolare ogni volta gli integrali.

La classe IntegraleMC:

Come al solito possiamo pensare ad una interfaccia generica IntegraleMC dalla quale derivano le classi concrete IntegratoreMedia e IntegratoreHoM

class IntegraleMC {

public:

  IntegraleMC(unsigned int seed) : 
    m_gen(seed)
  {
    m_errore = 0;
    m_punti  = 0;
  }

  virtual double Integra( const FunzioneBase& f, double inf , double sup , unsigned int punti, double fmax ) = 0 ;

  double GetErrore() const {return m_errore;}
  unsigned int GetN() const {return m_punti;}

protected:

  RandomGen m_gen;
  double m_errore;
  unsigned int m_punti;

};

class IntegratoreMedia : public IntegraleMC {

public:

  IntegratoreMedia(unsigned int seed) : IntegraleMC(seed) { ; };

  virtual double Integra ( const FunzioneBase& f, double inf , double sup , unsigned int punti, double fmax ) {

    // .... 

  };

};

notate che tra i data membri della classe compare un oggetto di tipo RandomGen m_gen: per inizializzare l'oggetto con il costruttore desiderato ( quello che accetta un unsigned int come seed) lo dobbiamo specificare nella lista di inizializzazione nella definizione del costruttore di IntegraleMC:

IntegraleMC(unsigned int seed) : m_gen(seed)  
  {

    // ....

  }

potete trovare una buona spiegazione dell'utilizzo delle liste di inizializzazione in qui.

Per salvare la struttura delle classi virtuale/concreta siamo stati costretti ad aggiungere il campo fmax anche ad IntegratoreMedia anche se non necessario. In questo caso non è veramente vantaggioso utilizzare questo tipo di schema.

Un vector di istogrammi:

In alcuni casi può essere comodo conservare una serie di istogrammi riempiti in un ciclo. Una possibilità è quella di creare puntatori a istogrammi e conservarli in un vector<TH1F*> per un eventuale uso seguente. Un esempio è tratteggiato qui sotto in caso sia utile :

  vector<TH1F*> vhistos;

  for ( int k = 0 ; k < cases.size() ; k++ ) {     // ciclo sui casi da studiare

    TH1F* h = new TH1F ( ... ) ;                   // costruzione istogramma 

    for ( int j = 0 ; j < n ; j++ ) h->Fill(...) ; // riempimento istogramma

    vhistos.push_back(h);                          // conserviamo i puntatori

  }

  for ( int k = 0 ; k < cases.size() ; k++ ) { 

    // ...

    vhistos[k]->Draw();

    // ...

  }

Calcolo di integrali con il metodo della media

Come discusso a lezione il metodo della media consiste nel valutare la media delle valutazioni della funzione $f(x)$ in un set di N punti $x_{i}$ casuali e distribuiti uniformemente tra $a$ (minimo estremo di integrazione) e $b$ (massimo estremo integrazione). La stima dell'integrale si ottiene poi dalla seguente formula: $\int_{a}^{b} f(x) dx \simeq (b-a) \frac{ \sum_{n=1}^{N} f(x_{n})}{N}$

Calcolo di integrali con il metodo hit-or-miss

Il metodo hit-or-miss si basa sulla generazione di una coppia di numeri $x \in [a,b]$ e $y \in [0,f_{max}]$ dove $f_{max}$ è un numero maggiore del massimo valore assunto da $f(x)$ nell'intervallo $[a,b]$ . Generata la coppia $x,y$ si incrementa un contatore $N_{tot}$ e si valuta quindi $f(x)$ : se $y < f(x)$ si incrementa anche il contatore $N_{hit}$ . La procedura viene ripetuta fino a che il numero di estrazioni è pari all' $N_{max}$ richiesto. La stima dell'integrale si ottiene poi dalla seguente formula: $\int_{a}^{b} f(x) dx \simeq (b-a) f_{max} \frac{N_{HIT}}{N_{TOT}}$

Risultati attesi

Le distribuzioni attese degli integrali per i diversi valori di N dovrebbero avere questo aspetto :

L'andamento dell'errore (deviazione standard delle distribuzioni precedenti) in funzione di N (in scala log-log):

La curva rossa rappresenta il risultato di un fit con una funzione del tipo $f(x) = k / \sqrt{N}$

Nel caso del metodo della media sono necessari ~ 900000 punti per ottenere un integrale con un errore di 0.001.

ESERCIZIO 10.3 - Calcolo di integrali multidimensionali con metodi MonteCarlo (facoltativo):

Provare a risolvere il seguente integrale utilizzando per esempio il metodo della media. $\int_{1}^{2} \int_{1}^{2} \left( \frac{5}{x^{2}} cos(4y) sin(7x) +10 \right) dxdy$ Suggerimento: si potrebbe costruire una classe FunzioneScalareBase astratta da cui la funzione integranda erediti. La classe `FunzioneScalareBase`` avrà un metodo virtuale puro

virtual double Eval(const vector<double>&) const=0;

Alla classe che calcola l'integrale si dovrà aggiungere un metodo dedicato del tipo

double Media(const FunzioneScalareBase* f, const vector<double>& inf, const vector& sup, unsigned int punti);

dove inf e sup sono i vettori che rappresentano gli estremi di integrazione.

Risultati attesi

Integrando la funzione con 100000 punti si dovrebbe ottenere un risultato $I = 10.241 +/- 0.004$

ESERCIZIO 10.4 - Errore nel caso di integrali multimensionali (facoltativo) :

Provare a ripetere le consegne dell'esercizio 10.2 applicate all'integrale multidimensionale dell'esercizio 10.3. In questo modo si può facilmente verificare che la legge con cui scala l'errore è indipendente dalla dimensione dell'integrale.

Qualche approfondimento su generatori di numeri casuali in C++11

Nel C++11 è stata inserita la libreria ufficiale per la generazione di numeri casuali: si veda per esempio qui. Provate a dare un'occhiata a questo codice per trovare un esempio su come utilizzare questa libreria e su come usare le librerie di ROOT ( si faccia riferimento alla referenza . )