Laboratorio Processi Stocastici Annalisa Pascarella Algoritmo istogramma INF

Laboratorio Processi Stocastici Annalisa Pascarella

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 VALORI = INF+DELTA/2 : DELTA : SUP-DELTA/2 figure bar(VALORI, contatore)

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

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)

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, 50) istogramma in 50 intervalli [counts bins] = hist(data, 50) i conteggi in counts, i punti medi degli intervalli in bins

Numeri casuali

Un po’ di storia n I numeri casuali sono utilizzati per costruire simulazioni di natura probabilistica di q q n 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

Un po’ di storia n 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 n Ulam fu uno dei personaggi chiave nel progetto americano per la costruzione della bomba atomica durante la II guerra mondiale q q 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

Cos’è un numero casuale? n n Lancio di un dado: l’imprevedibilità del numero ottenuto come punteggio conferisce allo stesso una forma di casualità Diversi metodi per generare numeri casuali q q 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 pseudocasuali ossia numeri generati da algoritmi numerici deterministici in grado di superare una serie di test statistici che conferiscono a tali numeri un’apparente casualità

Criteri n I fattori che determinano l’accettabilità di un metodo sono essenzialmente i seguenti: q q q 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.

Metodo middle-square n Genera numeri pseudo-casuali distribuiti in modo uniforme n In tale distribuzione uniforme ogni possibile numero in un determinato intervallo è ugualmente probabile q 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. Von Neumann)

Metodo middle-square n n n 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 q q q n eleviamo al quadrato e otteniamo le 8 (2 c) 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

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

Generatore lineare congruenziale Il metodo LCG ha bisogno di un seme per generare la n sequenza di numeri pseudo-casuiali secondo la seguente regola deterministica xn+1 q q n = (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 x 0 = 1 si ottiene per n=30 q 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

Bontà di un generatore LCG n Il problema della scelta dei migliori valori di a, c ed m è il punto cruciale del metodo q q n 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”: raggiungibile solo se 1. 2. 3. periodo 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 massimo M,

Bontà di un generatore LCG n Una delle scelte più popolari è q q 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

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

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

Verificare la casualità n n 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 n n 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 n n con scelta ottimale con i valori m=31, a=13, e c=0

Generatori e MATLAB n n 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 q q 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”

rand n n 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 n Per vedere gli algoritmi utilizzati da MATLAB q help rand

Esercizio n n 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 q q q 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.

Calcolo di p n n n 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

Esercizio n Calcolare p col metodo Monte Carlo q 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 q ripetere per diversi valori di N q q q

Aree e volumi n Il metodo Monte Carlo può essere usato anche per calcolare l’area della circonferenza n 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

Metodo Monte Carlo n n Vengono denominate le tecniche 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 U 1, …, Uk variabili casuali i. i. d. come U allora g(U 1), …, g(Uk) sono variabili casuali i. i. d. aventi come media q

Calcolo di integrali n 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) n Calcolare l’integrale q si prenda come X una v. a. normale standard e si usino n=10000

Calcolo di integrali n Calcolare l’integrale q n 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!)

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)

Metodo d’inversione. Il teorema ci fornisce una densità funzione di ripartizione 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.

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

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