Laboratorio Processi Stocastici

1. Laboratorio Processi Stocastici Annalisa Pascarella Istituto per le Applicazioni del Calcolo "M. Picone"Consiglio Nazionale delle Ricerche Roma

2. Informazioni • e-mail: • pascarella@dima.unige.it • annalisa.pascarella@unipr.it • a.pascarella@iac.cnr.it • webpage • www.dima.unige.it/~pascarel/ • Giovedì 10 Novembre PC2 12-14 • Venerdì 11 Novembre PC2 11-13 • Lunedì 5 Dicembre PC2 14-16

3. Programma • MATLAB • esercizi vari durante il laboratorio • MCMC • metodi Monte Carlo, • algoritmi per la generazione di numeri pseudo-casuali • Metropolis-Hastings • Simulazione processo di Poisson

4. MATLAB

5. MATLAB • MATrix LABoratory • Linguaggio di programmazione interpretato • legge un comando per volta eseguendolo immediatamente • Per avviarlo -> icona sul desktop workspace command window

6. MATLAB come calcolatrice 4 + 7 invio • è possibile definire variabili e operare su esse x = 9 -> invio Operatori aritmetici: + - * / ^ Caratteri speciali: ; % : Variabili predefinite: i, pi, NaN, Inf Funzioni elementari: sin, cos, log, exp help mean

7. Comandi utili • clear a • per cancellare una variabile dal workspace • clear all • per cancellare tutte le variabili dal workspace • ans • ultima variabile memorizzata • clc • pulisce lo schermo • help <nome_funzione>

8. Lavorare con MATLAB • In MATLAB tutte le variabili sono trattate come matrici • scalari -> matrici 1 x 1 • vettori riga -> matrici 1 x n v = (v1,…, vn) • vettori colonna -> matrici n x 1 v = (v1,…, vn)T • matrici -> matrici m x n

9. Vettori a = [1 2 3 4 5] a = [1, 2, 3, 4, 5] • Per definire un vettore riga • Per definire un vettore colonna • Usando : a = [1; 2; 3; 4; 5] a = [1 2 3 4 5] ’ a = 1:3:10 b = -5:5

10. Matrici A = [3 0; 1 2] A = [3 0 1 2] • Per definire una matrice • size(A) -> dimensioni della matrice • per memorizzare le dimensioni -> [r c] = size(A) b1 = [3;1] b2 = [0; 2] b3 = [3; 0] B = [b1, b2, b3] B(2,3) B(2,3) = 1; B • per selezionare un elemento • per modificare l’elemento • per visualizzare B

11. Il comando : • Importante per la manipolazione delle matrici • generazione di vettori che siano delle progressione aritmetiche di passo costante • a = [1:10] o a = 1:10 • b = 1: .2 : 4 • c = 3:0 -> non produce niente!!!! • c = 3: -1: 1 • mediante : si possono estrarre righe e colonne B(2,:) estrarre la riga R2 estrarre la colonna C2 B(:,2)

12. Identità-zero-uno eye(n) eye(3) identità di ordine n -> zeros(m,n) zeros(2,3) matrice nulla m x n -> • ones(m,n) • ones(2,3) matrice m x n di 1 ->

13. Operazioni size(A)= size(B) Somma / Differenza A+B, A-B A’ Trasposta #CA = #RB A*B Prodotto A*k Prodotto per uno scalare A.*B size(A) = size(B) Elemento per elemento

14. Script e funzioni • Script • parametri in ingresso non modificabili • le variabili usate sono messe nella memoria di lavoro di MATLAB • Funzioni • script al quale si possono passare parametri in ingresso ed ottenerne in uscita • sintassi • y1,…,yn-> parametri in uscita • x1,…,xn–> parametri in entrata • le variabili usate all’interno sono locali function [y1,…,yn] = nome_funzione(x1,…,xn)

15. Script • E’ possibile scrivere degli script in Matlab • cliccando su new • File -> New -> M-file

16. Le funzioni • L’m file va salvato col nome nome_funzione.m • il nome del file deve essere identico a quello della funzione • La funzione può essere richiamata • dalla finestra di comando • all’interno di uno script • da altre funzioni digitando [y1,…,yn]=nome_funzione(x1,…,xn) • Per poter richiamare la funzione dobbiamo essere nella directory nella quale è salvata la funzione oppure “settare” nel path di Matlab la directory nella quale la funzione è salvata.

17. Cicli for i = n1:passo:n2 blocco di istruzioni end • Ciclo incondizionato • Ciclo condizionato • Test condizionale while condizione blocco di istruzioni end if condizione1 blocco di istruzioni elseif condizione2 blocco di istruzioni else blocco di istruzioni end

18. Operatori • Operatori relazionali: • < , <= , > , >= , == , = , = • si usano per confrontare tra di loro gli elementi di 2 matrici; il risultato dell’operazione sarà • 0 se la relazione è falsa • 1 se la relazione è vera • Operatori logici: • & , | ,  • si usano per combinare tra loro gli operatori relazionali • Nota • = serve per assegnare valore ad una variabile • == per verificare se una variabile assume un determinato valore

19. Input\output • input • sprintf n = input(‘inserisci un intero ’); disp(sprintf(‘n = %d’,n)) disp(‘stringa di caratteri’)

20. Grafica • In MATLAB è possibile • disegnare funzioni in 2D e 3D • rappresentare graficamente dei dati • Il comando • si usa: • per rappresentare punti nel piano • per disegnare il grafico di una funzione • x e y devono essere vettori di ugual misura plot(x,y)

21. Esempio - I • Per rappresentare dei punti nel piano x = [1 2 3 7 -9 2]; y = [-2 -6 1 5 7 2]; plot(x,y) figure(2) plot(x,y,'*')

22. Esempio - II • Per “plottare” la funzione y=sin(x) definiamo l’intervallo in cui vogliamo disegnare la funzione x = [-pi:.01:pi]; y = sin(x); plot(x,y) definiamo la funzione disegniamo la funzione figure(2) plot(x,y, '-g') è possibile inserire un terzo parametro di input

23. Sintassi del comando “plot” plot(x, y) • x e y sono i vettori dei dati (ascisse e ordinate dei punti) • x e y come sopra; opzioni è una stringa opzionale che definisce il tipo di colore, di simbolo e di linea usato nel grafico. • help plot per vedere quali sono le varie opzioni • realizza il grafico del vettore y rispetto ai propri indici plot(x, y, 'opzioni') plot(y)

24. Comandi utili - I figure(num) • per creare (richiamare) una finestra grafica • per avere più grafici nella stessa finestra • holdoff per disattivare la funzione • per riscalare il grafico • per creare diversi grafici separati in una stessa finestra • esistono diversi comandi per “abbellire” i grafici • title, xlabel, ylabel, legend hold on axis([xminxmaxyminymax]) sublot(righe, colonne, sottofinestra)

25. Risultati usando hold on figure(1); hold on; grid on y2 = cos(x); plot(x,y2,’r’) title(‘seno e coseno’) % creiamo delle sottofinestre figure(3); subplot(1,2,1); plot(x,y); title('seno') subplot(1,2,2); plot(x,y2); title('coseno') usando subplot

26. Esercizio Scrivere una funzione che crei l’istogramma di un vettore Caricare il vettore dei dati nella variabile “data”: data = load(‘dato_per_istogramma.dat’); size(data) Osserviamo i dati plot(data) plot(data, ones(size(data)) , ’ . ’)

27. Algoritmo istogramma Scelta degli estremi e della larghezza intervallo INF DELTA SUP Contiamo quanti elementi del vettore cadono in ogni intervallo: creiamo un vettore il cui valore i-esimo rappresenti il numero di conteggi nell’i-esimo intervallo

28. Algoritmo istogramma • Per ogni elemento del vettore data(i) • Per ogni intervallo • Se data(i) è compreso nei valori dell’intervallo • Incrementare il contatore relativo a quell’intervallo INF DELTA SUP Il j-esimo intervallo ha come estremi INF+(j-1)*DELTA e INF+j*DELTA

29. 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)

30. 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

31. 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)

32. 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

33. Metodi Monte Carlo

34. 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, 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 • 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, uno dei personaggi chiave nel progetto americano per la costruzione della bomba atomica durante la II guerra mondiale)

35. 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à

36. 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. La generazione dei numeri casuali è troppo importante per essere lasciata al caso… (J.VonNeumann)

37. Generatori • Metodo middle-square • genera numeri pseudo-casuali distribuiti in modo uniforme • in tale distribuzione uniforme ogni possibile numero in un determinato intervallo è ugualmente probabile • 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

38. 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”

39. 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 • Una volta avviato MATLAB, il primo numero casuale generato è sempre lo stesso: 0.95012928514718 come anche la successione di numeri casuali • rand(‘state’,0)

40. 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

41. Metodo Monte Carlo • 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

42. Metodo Monte Carlo • L’idea è quella di estrarre un insieme i.i.d. di campioni da una pdf target p definita su uno spazio a grandi dimensioni ai quali sono associati dei pesi tale che l’integrale di una qualsiasi funzione misurabile rispetto alla pdf target p(x) possa essere approssimato dalla somma pesata • i pesi wisono determinati dalla stessa pdf • Tre approcci • randomsampling -> campiono direttamente dalla pdf target • importancesampling • MCMC

43. Random sampling • Se è un insieme i.i.d. di campioni generato dalla pdf target p il metodo Monte Carlo approssima la pdf target con la seguente funzione di densità empirica usando tale densità empirica si può calcolare un’approssimazione dell’integrale I per la legge dei grandi numeri si ha la convergenza a I(f)

44. Importancesampling • L’ipotesi principale per il randomsampling è che si sappia campionare da p(x), ma spesso si ha a che fare con pdf complicate. L’idea alla base dell’IS è usare una pdf dalla quale si sa campionare Se si possono estrarre N i.i.d. campioni da q(x) e calcolare i pesi p(x)/q(x) una stima Monte Carlo di I(f) sarà data da

45. Applicazioni • Viene utilizzato per: • simulazione di fenomeni naturali • simulazione di apparati sperimentali • calcolo di integrali • Problemi di natura statistica in cui Monte Carlo viene utilizzato per l’approssimazione di integrali sono ad esempio: • Inferenza Bayesiana → distribuzione a posteriori non appartiene a famiglie di distribuzioni note, dunque i momenti associati possono essere scritti sotto forma di integrale ma tipicamente non valutati analiticamente; • Problemi di massimizzazione della verosimiglianza → problemi inferenziali in cui la verosimiglianza stessa è funzione di uno o più integrali;

46. M/EEG MEG EEG Risoluzione temporale: 1 ms Campo magnetico [fT] Campo elettrico [mV]

47. Il problema inverso neuromagnetico • Il problema inverso della MEG/EEG consiste nel ricostruire l’evoluzione temporale delle sorgenti neuronali a partire dalle misure effettuate • Approccio statistico Bayesiano alla soluzione dei problemi inversi fisica matematici

48. probabilità a priori: contiene l’informazione che abbiamo a priori su J likelihood: contiene l’informazione sul problema diretto costante di normalizzazione Approccio Bayesiano • Ogni variabile è considerata come una variabile aleatoria (B e J sono le v.a. dei dati e dell’incognita mentre b e j sono le loro realizzazioni) • La soluzione del problema inverso è la densità di probabilità (pdf) dell’incognita condizionata dalle misure: Teoremadi Bayes::

49. Esercizio - 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

50. Esercizio - calcolo di p • 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