Title: Laboratorio Processi Stocastici
1 Laboratorio Processi Stocastici
- Annalisa Pascarella
- Istituto per le Applicazioni del Calcolo "M.
Picone"Consiglio Nazionale delle Ricerche - Roma
2Informazioni
- e-mail
- pascarella_at_dima.unige.it
- annalisa.pascarella_at_unipr.it
- a.pascarella_at_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
3Programma
- MATLAB
- esercizi vari durante il laboratorio
- MCMC
- metodi Monte Carlo,
- algoritmi per la generazione di numeri
pseudo-casuali - Metropolis-Hastings
- Simulazione processo di Poisson
4MATLAB
5MATLAB
- MATrix LABoratory
- Linguaggio di programmazione interpretato
- legge un comando per volta eseguendolo
immediatamente - Per avviarlo -gt
- icona sul desktop
workspace
command window
6MATLAB come calcolatrice
- è possibile definire variabili e operare su esse
Operatori aritmetici - / Caratteri
speciali
Variabili predefinite i, pi, NaN, Inf Funzioni
elementari sin, cos, log, exp
help mean
7Comandi 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 ltnome_funzionegt
8Lavorare con MATLAB
- In MATLAB tutte le variabili sono trattate come
matrici - scalari -gt matrici 1 x 1
- vettori riga -gt matrici 1 x n
- v (v1,, vn)
- vettori colonna -gt matrici n x 1
- v (v1,, vn)T
- matrici -gt matrici m x n
9Vettori
- 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 2 3 4 5 a 1 2 3 4 5
a 1310 b -55
10Matrici
A 3 0 1 2 A 3 0 1 2
- size(A) -gt dimensioni della matrice
- per memorizzare le dimensioni -gt r c
size(A)
b1 31 b2 0 2 b3 3 0 B b1, b2,
b3
- per selezionare un elemento
- per modificare lelemento
- per visualizzare B
11Il comando
- Importante per la manipolazione delle matrici
- generazione di vettori che siano delle
progressione aritmetiche di passo costante - a 110 o a 110
- b 1 .2 4
- c 30 -gt 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)
12Identità-zero-uno
eye(n) eye(3)
identità di ordine n -gt
zeros(m,n) zeros(2,3)
matrice nulla m x n -gt
ones(m,n) ones(2,3)
matrice m x n di 1 -gt
13Operazioni
size(A) size(B)
Somma / Differenza
AB, A-B
A
Trasposta
AB
CA RB
Prodotto
Ak
Prodotto per uno scalare
A.B
size(A) size(B)
Elemento per elemento
14Script 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 -gt parametri in uscita
- x1,,xn gt parametri in entrata
- le variabili usate allinterno sono locali
function y1,,yn nome_funzione(x1,,xn)
15Script
- E possibile scrivere degli script in Matlab
- cliccando su new
- File -gt New -gt M-file
16Le funzioni
- Lm 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
- allinterno di uno script
- da altre funzioni
- digitando y1,,ynnome_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.
17Cicli
for i n1passon2 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
18Operatori
- Operatori relazionali
- lt , lt , gt , gt , , , ?
- si usano per confrontare tra di loro gli elementi
di 2 matrici il risultato delloperazione 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
19Input\output
n input(inserisci un intero ) disp(sprintf(n
d,n))
disp(stringa di caratteri)
20Grafica
- 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)
21Esempio - 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,'')
22Esempio - II
- Per plottare la funzione ysin(x)
definiamo lintervallo in cui vogliamo disegnare
la funzione
x -pi.01pi 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
23Sintassi del comando plot
-
- 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)
plot(x, y, 'opzioni')
plot(y)
24Comandi utili - I
- per creare
(richiamare) una finestra grafica -
- per avere più
grafici nella stessa finestra - hold off 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
-
figure(num)
hold on
axis(xmin xmax ymin ymax)
sublot(righe, colonne, sottofinestra)
25Risultati
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
26Esercizio
Scrivere una funzione che crei listogramma di un
vettore
Caricare il vettore dei dati nella variabile
data data load(dato_per_istogramma.dat) s
ize(data)
Osserviamo i dati plot(data) plot(data,
ones(size(data)) , . )
27Algoritmo istogramma
Scelta degli estremi e della larghezza intervallo
INF
SUP
DELTA
Contiamo quanti elementi del vettore cadono in
ogni intervallo creiamo un vettore il cui valore
i-esimo rappresenti il numero di conteggi
nelli-esimo intervallo
28Algoritmo istogramma
Per ogni elemento del vettore data(i) Per ogni
intervallo Se data(i) è compreso nei valori
dellintervallo Incrementare il contatore
relativo a quellintervallo
INF
SUP
DELTA
Il j-esimo intervallo ha come estremi
INF(j-1)DELTA e INFjDELTA
29Algoritmo 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 1size(data,2) per ogni
dato for j 1 NUM_INT per ogni intervallo
if data(i)gtINF(j-1)DELTA
data(i)ltINFjDELTA contatore(j)
contatore(j)1 end end end VALORI INFDELTA/2
DELTA SUP-DELTA/2 figure bar(VALORI,
contatore)
30Algoritmo istogramma (efficiente)
Lalgoritmo appena scritto fa un ciclo di
troppo...
INF
SUP
1
2
k
Osserviamo che il singolo valore data(i) INF lt
data(i) lt SUP 0 lt data(i)-INF lt
SUP-INFDELTANUM_INT 0 lt (data(i)-INF)/DELTA lt
NUM_INT
31Algoritmo 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 1size(data,2) per ogni
dato j ceil((data(i)-INF)/DELTA) contatore(j)
contatore(j) 1 end VALORI INFDELTA/2
DELTA SUP-DELTA/2 figure bar(VALORI, contatore)
32Istogrammi e MATLAB
Esiste un comando che fa listogramma 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
33Metodi Monte Carlo
34Un 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 - lidea 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)
35Cosè un numero casuale?
- Lancio di un dado limprevedibilità 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 unapparente casualità
36Criteri
- I fattori che determinano laccettabilità 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
dallelaboratore e deve consumare poco spazio di
memoria.
La generazione dei numeri casuali è troppo
importante per essere lasciata al caso (J.Von
Neumann)
37Generatori
- 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 - xn1 (axnc)mod m , ngt0
- con a,c ed m opportuni numeri interi costanti
- xn1 assume valori compresi tra 0, , m-1
38Generatori 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 dellesponente con la data
della scoperta dellAmerica questo generatore è
comunemente chiamato il generatore di Cristoforo
Colombo
39rand
- La funzione rand genera una successione di numeri
casuali distribuiti uniformemente nellintervallo
(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)
40Metodo Monte Carlo
- Vengono denominate le tecniche che utilizzano
variabili casuali per risolvere vari problemi,
anche non di natura aleatoria. - Vediamo lapproccio generale supponiamo che un
problema si riconduca al calcolo di un integrale -
- Sia U la variabile casuale uniforme, allora
41Metodo 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 -
42Metodo Monte Carlo
- Lidea è 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 lintegrale di
una qualsiasi funzione misurabile rispetto alla
pdf target p(x) - possa essere approssimato dalla somma pesata
- i pesi wi sono determinati dalla stessa pdf
- Tre approcci
- random sampling -gt campiono direttamente dalla
pdf target - importance sampling
- MCMC
43Random 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 unapprossimazione dellintegrale I - per la legge dei grandi numeri si ha la
convergenza a I(f)
44Importance sampling
- Lipotesi principale per il random sampling è che
si sappia campionare da p(x), ma spesso si ha a
che fare con pdf complicate. Lidea alla base
dellIS è 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 -
45Applicazioni
- 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 lapprossimazione 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
46M/EEG
MEG
EEG
Risoluzione temporale 1 ms
Campo magnetico fT
Campo elettrico mV
47Il problema inverso neuromagnetico
fisica
matematici
- Il problema inverso della MEG/EEG consiste nel
ricostruire levoluzione temporale delle sorgenti
neuronali a partire dalle misure effettuate - Approccio statistico Bayesiano alla soluzione dei
problemi inversi
48Approccio Bayesiano
- Ogni variabile è considerata come una variabile
aleatoria (B e J sono le v.a. dei dati e
dellincognita mentre b e j sono le loro
realizzazioni) - La soluzione del problema inverso è la densità di
probabilità (pdf) dellincognita condizionata
dalle misure
Teorema di Bayes
49Esercizio - 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 allinterno 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 larea della circonferenza e
- quella del quadrato
- può essere usato per stimare p
50Esercizio - calcolo di p
- Calcolare p col metodo Monte Carlo
- considerare un quadrato di lato 2 (come in
figura) il cui centro coincide con lorigine 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 allinterno del cerchio - stimare p usando la formula
-
- ripetere per diversi valori di N
51Campionamento
- In generale non è sufficiente utilizzare sequenze
di numeri casuali distribuiti uniformemente - in molti problemi è necessario disporre di numeri
casuali estratti da densità di probabilità
diverse, quali la normale, lesponenziale, la
poissoniana, etc - Esempio
- Simulare lenergia di una particella con una
distribuzione gaussiana intorno ad un valore
medio e con una data sigma - Varie possibilità
- Metodo dellinversione
- Metodo del rigetto
52Generare 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
nellintervallo 0,1 Quindi per campionare una
variabile aleatoria X con distribuzione F basta
campionare una variabile uniforme in 0,1 e poi
considerare XF-1(U)
53Metodo dinversione
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
54Esempio 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
55In MALTAB...
Ora provate... data rand(1,1000) hist(data) da
ta exprand(1,1,1000) hist(data) poissrnd Poiss
on randn Gaussiana
56Variabile aleatoria discreta
- Supponiamo di voler simulare una variabile
aleatoria discreta X che può assumere m valori
distinti xi, i1,,m con probabilità - Sia Fk la probabilità cumulativa
- Si verifichi che se U è una v.a. uniforme in
0,1 allora la variabile - ha la probabilità desiderata
57Variabile aleatoria discreta
- Si scriva una funzione Matlab che, dati un intero
n gt 0 e i vettori x x1, x2, . . . , xm e p
p1, p2, . . . , pm, restituisce il vettore y
contenente n campionamenti della variabile
aleatoria discreta X. - Si consideri la variabile aleatoria X tale che
- si calcolino 1000 campioni della variabile
aleatoria e si verifichi la correttezza dei
risultati ottenuti confrontando media, varianza e
funzione di ripartizione cumulativa (cdf) con i
valori teorici.
58Metodo del rigetto
- Obiettivo Estrarre una sequenza di numeri
casuali secondo la distribuzione f(x),
nellintervallo (x1, x2) - Metodo
- Generare x uniforme in (x1,x2)
- Generare y uniforme in (0,fmax)
- Valutare f(x)
- Confrontare y con f(x)
- Se y lt f(x) accettare x
- Se y gt f(x) ripetere da a) in poi
fmax
x1
x2
59Metodo del rigetto
- Si può migliorare lefficienza del metodo
effettuando il campionamento non in un
rettangolo ma in una regione definita da una
funzione g(x) maggiorante di f(x). - Se si è in grado di generare numeri casuali
distribuiti secondo g(x) per esempio con il
metodo dellinversione, la procedura
dellinversione diventa - Si estrae x secondo g(x)
- Si estrae u con distribuzione uniforme
- tra 0 e g(x)
- Se ultf(x) si accetta x
60Metodo del rigetto
- I presupposti per utilizzare questo metodo sono
costruire unopportuna nuova distribuzione di
probabilità g da cui sappiamo come campionare, e
definire una funzione envelope e(x) tale che
e(x) g(x)/a f (x) - Campioniamo Y dalla distribuzione g
- Campioniamo U Unif (0, 1)
- Eliminiamo il valore Y se U lt f (Y)/e(Y)
- Altrimenti, manteniamo il valore XY come un
campionamento dalla distribuzione obiettivo f - torniamo al punto iniziale.
- Se riprendiamo il terzo punto, è come se avessimo
campionato da Unif (0, e(y)) ed accettassimo il
valore se risulta essere inferiore ad f (y)