Laboratorio Processi Stocastici - PowerPoint PPT Presentation

1 / 60
About This Presentation
Title:

Laboratorio Processi Stocastici

Description:

Programma. MATLAB. esercizi vari durante il laboratorio. MCMC. metodi Monte Carlo, algoritmi per la generazione di numeri pseudo-casuali. Metropolis-Hastings – PowerPoint PPT presentation

Number of Views:106
Avg rating:3.0/5.0
Slides: 61
Provided by: unig155
Category:

less

Transcript and Presenter's Notes

Title: 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_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

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 -gt
  • icona sul desktop

workspace
command window
6
MATLAB 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
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 ltnome_funzionegt

8
Lavorare 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

9
Vettori
  • 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
10
Matrici
A 3 0 1 2 A 3 0 1 2
  • Per definire una matrice
  • 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

11
Il 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)
12
Identità-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
13
Operazioni
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
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 -gt parametri in uscita
  • x1,,xn gt parametri in entrata
  • le variabili usate allinterno sono locali

function y1,,yn nome_funzione(x1,,xn)
15
Script
  • E possibile scrivere degli script in Matlab
  • cliccando su new
  • File -gt New -gt M-file

16
Le 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.

17
Cicli
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
18
Operatori
  • 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

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 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
23
Sintassi 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)
24
Comandi 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)
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 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)) , . )
27
Algoritmo 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
28
Algoritmo 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
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 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)
30
Algoritmo 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
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 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)
32
Istogrammi 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
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
  • 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)

35
Cosè 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à

36
Criteri
  • 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)
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
  • xn1 (axnc)mod m , ngt0
  • con a,c ed m opportuni numeri interi costanti
  • xn1 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 dellesponente con la data
    della scoperta dellAmerica questo generatore è
    comunemente chiamato il generatore di Cristoforo
    Colombo

39
rand
  • 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)

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

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
  • 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

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 unapprossimazione dellintegrale I
  • per la legge dei grandi numeri si ha la
    convergenza a I(f)

44
Importance 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

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

46
M/EEG
MEG
EEG
Risoluzione temporale 1 ms
Campo magnetico fT
Campo elettrico mV
47
Il 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

48
Approccio 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
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 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

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

51
Campionamento
  • 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

52
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
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)
53
Metodo 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
54
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
55
In MALTAB...
Ora provate... data rand(1,1000) hist(data) da
ta exprand(1,1,1000) hist(data) poissrnd Poiss
on randn Gaussiana
56
Variabile 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

57
Variabile 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.

58
Metodo 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
59
Metodo 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

60
Metodo 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)
Write a Comment
User Comments (0)
About PowerShow.com