Laboratorio Processi Stocastici

1. Laboratorio Processi Stocastici Annalisa Pascarella

2. Algoritmo istogramma • INF = -4; • SUP = 4; • DELTA = 0.4; • NUM_INT = (SUP-INF)/DELTA; % numero di intervalli • contatore = zeros(1,NUM_INT) % inizializziamo il contatore; • for i = 1:size(data,2) % per ogni dato • for j = 1: NUM_INT % per ogni intervallo • if data(i)>INF+(j-1)*DELTA && data(i)<INF+j*DELTA • contatore(j) = contatore(j)+1; • end • end • end • VALORI = INF+DELTA/2 : DELTA : SUP-DELTA/2 • figure • bar(VALORI, contatore)

3. Algoritmo istogramma (efficiente) L’algoritmo appena scritto fa un ciclo di troppo... INF SUP 1 2 k Osserviamo che il singolo valore data(i) INF < data(i) < SUP 0 < data(i)-INF < SUP-INF=DELTA*NUM_INT 0 < (data(i)-INF)/DELTA < NUM_INT

4. Algoritmo istogramma (efficiente) • INF = -4; • SUP = 4; • DELTA = 0.4; • NUM_INT = (SUP-INF)/DELTA; % numero di intervalli • contatore = zeros(1,NUM_INT) % inizializziamo il contatore; • for i = 1:size(data,2) % per ogni dato • j = ceil((data(i)-INF)/DELTA); • contatore(j) = contatore(j) + 1; • end • VALORI = INF+DELTA/2 : DELTA : SUP-DELTA/2 • figure • bar(VALORI, contatore)

5. Istogrammi e MATLAB Esiste un comando che fa l’istogramma delle frequenze dei valori di un vettore data = load(‘dato_per_istogramma.dat’) hist(data) hist(data,50) istogramma in 50 intervalli [counts bins] = hist(data,50) i conteggi in counts, i punti medi degli intervalli in bins

6. Numeri casuali

7. Un po’ di storia • I numeri casuali sono utilizzati per costruire simulazioni di natura probabilistica di • fenomeni fisici: reattori nucleari, traffico stradale, aerodinamica • problemi decisionali e finanziari: econometria, previsione Dow-Jones • informatica: rendering • varia natura: videogiochi • Il legame che esiste tra il gioco e le simulazioni probabilistiche è sottolineato dal fatto che a tali simulazioni è dato il nome di metodi Monte Carlo

8. Un po’ di storia • L’idea di utilizzare in modo sistematico simulazioni di tipo probabilistico per risolvere un problema di natura fisica viene generalmente attribuita al matematico polacco Ulam • Ulam fu uno dei personaggi chiave nel progetto americano per la costruzione della bomba atomica durante la II guerra mondiale • il progetto richiedeva la risoluzione di un enorme numero di problemi incredibilmente complessi • l’idea di utilizzare simulazioni casuali per risolvere tali problemi gli venne giocando a carte

9. Cos’è un numero casuale? • Lancio di un dado: l’imprevedibilità del numero ottenuto come punteggio conferisce allo stesso una forma di casualità • Diversi metodi per generare numeri casuali • hardware • calcolatore: il calcolatore è un oggetto puramente deterministico e quindi prevedibile, per cui nessun calcolatore è in grado di generare numeri puramente casuali, ma solo numeri pseudo-casuali ossia numeri generati da algoritmi numerici deterministici in grado di superare una serie di test statistici che conferiscono a tali numeri un’apparente casualità

10. Criteri • I fattori che determinano l’accettabilità di un metodo sono essenzialmente i seguenti: • i numeri della sequenza generata devono essere uniformemente distribuiti (cioè devono avere la stessa probabilità di presentarsi); • i numeri devono risultare statisticamente indipendenti; • la sequenza deve poter essere riprodotta; • la sequenza deve poter avere un periodo di lunghezza arbitraria; • il metodo deve poter essere eseguito rapidamente dall’elaboratore e deve consumare poco spazio di memoria.

11. Metodo middle-square • Genera numeri pseudo-casuali distribuiti in modo uniforme • In tale distribuzione uniforme ogni possibile numero in un determinato intervallo è ugualmente probabile • ad es. se lanciamo un dato un certo numero di volte ognuna delle facce da 1 a 6 si presenterà circa 1/6 delle volte originando così una successione uniforme di numeri casuali compresi tra 1 e 6 La generazione dei numeri casuali è troppo importante per essere lasciata al caso… (J.VonNeumann)

12. Metodo middle-square • Supponiamo di voler generare un numero casuale di 4 cifre • Il metodo richiede come tutti i generatori di numeri casuali un valore iniziale, detto seme, dal quale vengono generati i successivi valori • Ad es. a partire da 1234 avente 4(c) cifre • eleviamo al quadrato e otteniamo le 8 (2c) cifre 01522756 • consideriamo solo le 4 (c) cifre di mezzo 5227 • ripetiamo il procedimento ottenendo 27321529 e 3215 e così via • Ogni nuovo numero è determinato univocamente dal predecessore. Ogni successione di numeri generata da questo algoritmo si ripeterà prima o poi. Il numero di numeri della sequenza prima che intervenga una ripetizione è detta periodo della sequenza

13. Esempio • Simulazione del lancio di un dado • definiamo il risultato ottenuto come d =1+[5ms/10^4] dove ms è il numero generato tramite il metodo middle-square • simulando 10 lanci consecutivi a partire dal seme 8022 otteniamo risultati che sembrano abbastanza realistici 5 3 3 4 3 3 4 2 3 1 • basta aumentare il numero di lanci per ottenere risultati non soddisfacenti (la successione ha periodo 38)

14. Generatore lineare congruenziale • Il metodo LCG ha bisogno di un seme per generare la sequenza di numeri pseudo-casuiali secondo la seguente regola deterministica xn+1 = (axn+c)mod m , n>=0 • con a,c ed m opportuni numeri interi costanti • xn+1 assume valori compresi tra 0, …, m-1 • Ad es. per a=13, c=0 (generatore puramente moltiplicativo) ed m=31 partendo da x0 = 1 si ottiene per n=30 1 13 14 27 10 6 16 22 7 29 5 3 8 11 19 30 18 17 4 21 25 15 9 24 2 26 28 23 20 12 • tale successione ha periodo 30 (= m-1). tutti i numeri da 1 a 30 compaiono per poi ripetersi

15. Bontà di un generatore LCG • Il problema della scelta dei migliori valori di a, c ed m è il punto cruciale del metodo • un aspetto importante è la lunghezza del periodo che dovrà essere molto grande, per cui m dovrà essere grande • un altro aspetto consiste nel garantire che per un dato m i valori di a, c siano tali che la successione abbia periodo massimo • Generatore “periodico”: periodo massimo M, raggiungibile solo se • c e M sono primi tra loro • a-1 è divisibile per tutti i fattori primi di M • a-1 è multiplo di 4 se M è multiplo di 4

16. Bontà di un generatore LCG • Una delle scelte più popolari è m=231-1, a=75, c=0 • questo garantisce un periodo di 231-2=2147483646 ossia oltre 2 miliardi di numeri pseudo-casuali • il fatto che 231-1 sia un numero primo è fondamentale al fine di ottenere il massimo periodo xn+1 = (axn+c)mod m , n>=0

17. Un algoritmo per generare numeri random a = 7^5 M = 2^(31)-1 c=0 L(1) = 1; for i = 2:100 L(i) = mod(a*L(i-1)+c , M) u(i) = L(i)/M end resto = mod(dividendo,divisore) Gli u(i) sono distribuiti in maniera uniforme tra 0 e 1. Provare per credere

18. Verifica funzionamento 1. Fare istogramma dei numeri random generati 2. Modificare la lunghezza del vettore di numeri casuali (ad es. 100, 1,000 e 10,000) e osservare la “omogeneità” della distribuzione

19. Verificare la casualità • Una richiesta importante al fine di valutare la bontà di un generatore uniforme di numeri pseudo-casuali è l’assenza di correlazione tra i numeri generati dell’algoritmo. Non deve emergere nessuna relazione tra xn e xn+1 per n>0. • Questa proprietà può essere verificata graficamente realizzando il grafico di (xn, xn+j) per j>0 • nel grafico non dovranno comparire linee, forme o altre strutture regolari • Provare a disegnare il grafico per j=1 con 1000 punti ottenuti con il generatore LCG • con scelta ottimale • con i valori m=31, a=13, e c=0

20. Generatori e MATLAB • I generatori di numeri casuali più recenti non sono basati sul metodo LCG, ma sono una combinazione di operazioni di spostamento di registri e manipolazione sui bit che non richiedono nessuna operazione di moltiplicazione o divisione. Questo nuovo approccio risulta estremamente veloce e garantisce periodi incredibilmente lunghi • Nelle ultime versioni di MATLAB il periodo è 21492 • un milione di numeri casuali al secondo richiederebbe 10435 anni prima di ripetersi! • data la coincidenza dell’esponente con la data della scoperta dell’America questo generatore è comunemente chiamato il “generatore di Cristoforo Colombo”

21. rand • La funzione rand genera una successione di numeri casuali distribuiti uniformemente nell’intervallo (0,1) • La sintassi di tale funzione è rand(n,m) che genera una matrice n x m di numeri casuali distribuiti uniformemente • Per vedere gli algoritmi utilizzati da MATLAB • help rand

22. Esercizio • Sia X una v.a. uniforme nell’intervallo [0,1]. La si campioni n volte, con n=102, 103, 104, 105. • Per ciascun valore di n • si calcolino media e varianza campionarie (mediante i comandi mean e var) • si visualizzi l’istogramma dei valori campionati • si visualizzi la funzione di ripartizione empirica dei dati mediante il comando cdfplot e la si confronti graficamente con la funzione di ripartizione cumulativa di X.

23. Calcolo di p • Supponiamo di lanciare N freccette ad un bersaglio formato da un quadrato di lato L contenente una circonferenza • Assumiamo che le freccette siano lanciate casualmente all’interno del quadrato e che quindi colpiscano il quadrato in ogni posizione con uguale probabilità • Dopo molti lanci la frazione di freccette che ha colpito la circonferenza sarà uguale al rapporto tra l’area della circonferenza e quella del quadrato può essere usato per stimare p

24. Esercizio • Calcolare p col metodo Monte Carlo • considerare un quadrato di lato 2 (come in figura) il cui centro coincide con l’origine di un sistema di riferimento Oxy e una circonferenza inscritta in esso • generare 2 vettori, x e y, di numeri casuali di lunghezza N • calcolare il numero dei punti (NC) (x,y) così generati che cadono all’interno del cerchio • stimare p usando la formula • ripetere per diversi valori di N

25. Aree e volumi • Il metodo Monte Carlo può essere usato anche per calcolare l’area della circonferenza • La generalizzazione al calcolo di volumi nello spazio è immediata. Indicato con L il lato del cubo contenente la figura di cui si vuole misurare il volume V avremo dove Nc è il numero di punti generati in modo uniforme nel cubo e interni alla figura di cui si vuole misurare il volume

26. Metodo Monte Carlo • Vengono denominate le tecniche che utilizzano variabili casuali per risolvere vari problemi, anche non di natura aleatoria. • Vediamo l’approccio generale: supponiamo che un problema si riconduca al calcolo di un integrale Sia U la variabile casuale uniforme, allora Siano U1, …, Uk variabili casuali i.i.d. come U allora g(U1), …, g(Uk) sono variabili casuali i.i.d. aventi come media q

27. Calcolo di integrali • Sia X la v.c. avente densità p e Y la v.a. Y=f(X). Il valore Supponendo di essere capaci di campionare X l’integrale I può essere approssimato mediante il metodo Monte Carlo utilizzando lo stimatore di media campionaria per la v.a. Y=f(X) • Calcolare l’integrale • si prenda come X una v.a. normale standard e si usino n=10000 campioni

28. Calcolo di integrali • Calcolare l’integrale • si prenda come X una v.a. normale standard e si usino n=10000 campioni • Il metodo fornisce la stima 8.6080 con un errore del 4% circa (avendo usato ben 104 campioni!)

29. Generare numeri casuali con distribuzione arbitraria • Metodo di inversione • Sia X una variabile aleatoria continua a valori in R e F: (0,1)  R , la corrispondente funzione di ripartizione cumulativa: • La variabile aleatoria U = F(X) ha una densità di probabilità uniforme nell’intervallo [0,1] • Quindi per campionare una variabile aleatoria X con distribuzione F basta campionare una variabile uniforme in [0,1] e poi considerare X=F-1(U)

30. Metodo d’inversione Il teorema ci fornisce una regola per generare numeri con distribuzione arbitraria: se conosciamo F, prendiamo i numeri {ui} distribuiti secondo la legge uniforme e {F-1(ui)} sono distribuiti secondo F. densità funzione di ripartizione

31. Esempio: distribuzione esponenziale Generare numeri distribuiti secondo la legge esponenziale: se i numeri {ui} sono distribuiti secondo la legge uniforme, {F-1(ui)} hanno F come funzione di ripartizione. La variabile X~exp(l) ha funzione di ripartizione • La variabile X può essere ottenuta come trasformazione di una variabile uniforme

32. In MALTAB... Ora provate... data = rand(1,1000) hist(data) data = exprand(1,1,1000) hist(data) poissrnd Poisson randn Gaussiana