laboratorio processi stocastici n.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
Laboratorio Processi Stocastici PowerPoint Presentation
Download Presentation
Laboratorio Processi Stocastici

Loading in 2 Seconds...

play fullscreen
1 / 60

Laboratorio Processi Stocastici - PowerPoint PPT Presentation


  • 115 Views
  • Uploaded on

Laboratorio Processi Stocastici. Annalisa Pascarella Istituto per le Applicazioni del Calcolo "M. Picone " Consiglio Nazionale delle Ricerche Roma. Informazioni. e-mail: pascarella@dima.unige.it annalisa.pascarella@unipr.it a.pascarella@iac.cnr.it webpage

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'Laboratorio Processi Stocastici' - parson


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
laboratorio processi stocastici

Laboratorio Processi Stocastici

Annalisa Pascarella

Istituto per le Applicazioni del Calcolo "M. Picone"Consiglio Nazionale delle Ricerche

Roma

informazioni
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
programma
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
matlab1
MATLAB
  • MATrix LABoratory
  • Linguaggio di programmazione interpretato
    • legge un comando per volta eseguendolo immediatamente
  • Per avviarlo ->

icona sul desktop

workspace

command

window

matlab come calcolatrice
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

comandi utili
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>
lavorare con matlab
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
vettori
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

matrici
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
il comando
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)

identit zero uno
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 ->

operazioni
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

script e funzioni
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)

script
Script
  • E’ possibile scrivere degli script in Matlab
    • cliccando su new
    • File -> New -> M-file
le funzioni
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.
cicli
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

operatori
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
input output
Input\output
  • input
  • sprintf

n = input(‘inserisci un intero ’);

disp(sprintf(‘n = %d’,n))

disp(‘stringa di caratteri’)

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

esempio i
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,'*')

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

sintassi del comando plot
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)

comandi utili i
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)

risultati
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

esercizio
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)) , ’ . ’)

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

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

algoritmo istogramma2
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)
algoritmo istogramma efficiente
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 efficiente1
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
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

un po di storia
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)
cos un numero casuale
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à
criteri
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)

generatori
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
generatori e matlab
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”
slide39
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)
metodo monte carlo
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

metodo monte carlo1
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
metodo monte carlo2
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
random sampling
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)

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

applicazioni
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;
m eeg
M/EEG

MEG

EEG

Risoluzione temporale: 1 ms

Campo magnetico [fT]

Campo elettrico [mV]

il problema inverso neuromagnetico
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

approccio bayesiano

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

esercizio calcolo di p
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

esercizio calcolo di p1
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
campionamento
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, l’esponenziale, la poissoniana, etc
  • Esempio:
    • Simulare l’energia di una particella con una distribuzione gaussiana intorno ad un valore medio e con una data sigma
  • Varie possibilità:
    • Metodo dell’inversione
    • Metodo del rigetto
generare numeri casuali con distribuzione arbitraria
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 considerareX=F-1(U)
metodo d inversione
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

esempio distribuzione esponenziale
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
In MALTAB...

Ora provate...

data = rand(1,1000)

hist(data)

data = exprand(1,1,1000)

hist(data)

poissrnd Poisson

randn Gaussiana

variabile aleatoria discreta
Variabile aleatoria discreta
  • Supponiamo di voler simulare una variabile aleatoria discreta X che può assumere m valori distinti xi, i=1,…,m con probabilità

Sia Fkla probabilità cumulativa

Si verifichi che se U è una v.a. uniforme in [0,1] allora la variabile

ha la probabilità desiderata

variabile aleatoria discreta1
Variabile aleatoria discreta
  • Si scriva una funzione Matlab che, dati un intero n > 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.

metodo del rigetto
Metodo del rigetto
  • Obiettivo: Estrarre una sequenza di numeri casuali secondo la distribuzione f(x), nell’intervallo (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 < f(x) accettare x
      • Se y > f(x) ripetere da a) in poi

fmax

x1

x2

metodo del rigetto1
Metodo del rigetto
  • Si può migliorare l’efficienza 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 dell’inversione, la procedura dell’inversione diventa:
    • Si estrae x secondo g(x)
    • Si estrae u con distribuzione uniforme

tra 0 e g(x)

    • Se u<f(x) si accetta x
metodo del rigetto2
Metodo del rigetto
  • I presupposti per utilizzare questo metodo sono: costruire un’opportuna nuova distribuzione di probabilità g da cui sappiamo come campionare, e definire una funzione “envelope” e(x) tale che e(x) = g(x)/α ≥ f (x)
    • Campioniamo Y dalla distribuzione g
    • Campioniamo U ∼ Unif (0, 1)
    • Eliminiamo il valore Y se U < f (Y)/e(Y)
    • Altrimenti, manteniamo il valore X=Y 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)