Diapositiva 1 - Istituto per le Applicazioni del Calcolo

Laboratorio
Processi Stocastici
Annalisa Pascarella
Istituto per le Applicazioni del Calcolo "M. Picone"
Consiglio Nazionale delle Ricerche
Metodo Monte Carlo

N
L’idea è quella di estrarre un insieme i.i.d. di campioni X i i 1 da
una pdf target p definita su uno spazio a grandi dimensioni ai
quali sono associati dei pesi wi iN1 tale che l’integrale di una
qualsiasi funzione misurabile rispetto alla pdf target p(x)
E[ f ( X )]  I ( f ) 
X f ( x) p( x)dx
possa essere approssimato dalla somma pesata
N


I N ( f )   wi f ( X i )
i 1
i pesi wi sono determinati dalla stessa pdf
Tre approcci



random sampling -> campiono direttamente dalla pdf target
importance sampling -> uso una pdf dalla quale so campionare
MCMC
Random sampling

Se X i iN1 è 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
1 N
p N ( x )    ( x  xi )
N i 1
Usando
tale
densità
empirica
un’approssimazione dell’integrale I
si
può
1 N
I N ( f )   f ( xi )
N i 1

calcolare
Per la legge dei grandi numeri si ha la convergenza a I(f)
La velocità di convergenza dipende da N e non dallo spazio nel
quale vivono i campioni

principale vantaggio rispetto alle tecniche di quadratura
Random sampling

I campioni così ottenuti possono essere usati per calcolare ad es
xˆ  arg max p ( xi )

Se ad es f(x)=x e la pdf target è la pdf a posteriori p(x|y), I(f) è il
valor medio condizionale e una sua stima è
xCM

1 N
  xi
N i 1
La principale hp è che si possano estrarre N i.i.d. campioni dalla
pdf target


nel caso gaussiano esistono diverse routine per campionare da esse
in generale si devono campionare pdf complicate => per ottenere un
insieme di campioni con cui approssimare l’integrale I(f) si ricorre a
tecniche più sofisticate come l’importance sampling o i metoti MCMC
Importance sampling

L’idea alla base dell’IS è usare una pdf q(x) dalla quale si sa
campionare
I( f ) 




X
f ( x)
p( x)
p( x)
q ( x)dx  E q { f ( x)
}  E q { f ( x) w( x)}
q( x)
q( x)
q(x) prende il nome di proposal pdf
w(x) = p(x)/q(x)
il supporto di q deve contenere quello della pdf target p
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

1 N
I N ( f )   f ( xi )w( xi )
N i 1
in pratica stiamo effettuando un random sampling di f(x)w(x)
Importance sampling

Se sono in un contesto di inferenza Bayesiana vorrei poter
campionare p(x|y)


p ( y | x)p pr ( x)
p post ( x | y ) 
p ( y)
campionare tale pdf usando l’IS non è sempre possibile in quanto spesso
non si può calcolare la costante di normalizzazione
Una possibile soluzione


Markox Chain Monte Carlo: usano catene di Markov per generare
campioni da usare nell’integrazione Monte Carlo.
L’algoritmo Metropolis è considerato tra uno dei dieci metodi che hanno
avuto maggiore influenza sulle scienze e l’ingegneria nel XX secolo
MCMC
Da MC a MCMC

L’obiettivo dei metodi MC, e il loro principale utilizzo nelle
applicazioni, è l’integrazione



stima del valore atteso di una funzione di X ∼ p (x) tramite simulazione di
un campione da p (distribuzione target che può essere valutata facilmente
ma dalla quale non è immediato campionare).
simulare un campione da una generica funzione target p può risultare di
fondamentale importanza anche in altri contesti
Il metodo Markov Chain Monte Carlo (MCMC), che ci
consente di ottenere campioni da funzioni target qualsiasi,
consente quindi di raggiungere obiettivi inferenziali più generali
del solo calcolo di un integrale.
MCMC


L’idea di base dei metodi MCMC è di generare un campione
dalla distribuzione d’interesse p costruendo una catena di
Markov irriducibile e aperiodica, avente la distribuzione target
p come distribuzione stazionaria; per n sufficientemente
grande, una realizzazione Xn della catena è equivalente ad un
campionamento da p .
L’applicazione più popolare dei metodi MCMC è nell’ambito
dell’inferenza Bayesiana, dove la distribuzione target p è la
distribuzione a posteriori, generalmente indicata con π, ed X
sono i parametri di interesse, generalmente indicati con θ.
Differenze

La differenza sostanziale tra il metodo Monte Carlo classico e i
metodi MCMC è che nel primo caso i campioni generati sono
indipendenti tra loro, mentre nel secondo caso vengono generati
attraverso un processo stocastico di Markov.


nei metodi MCMC i campioni sono correlati tra loro e di conseguenza vi è
la necessità di un maggior numero di iterazioni per avere un risultato
sufficientemente accurato.
Non sempre è possibile trovare una distribuzione da cui generare
i campioni indipendenti, mentre è molto semplice generare una
catena di Markov che si muova nello spazio delle soluzioni,
addirittura anche nel caso in cui le probabilità che stiamo
cercando siano proporzionali a una costante a noi sconosciuta, o
di cui è difficile trovare il valore esatto.
Catene di Markov

Sia Xt il valore di una v.a. a t e sia S={s1, …, sm} l’insieme degli
stati. Il processo stocastico {Xt} è un processo di Markov se
Pr( X t 1  s j | X 0  s k ,..., X t  si )  Pr( X t 1  s j | X t  si )
Catene di Markov


Nella dinamica delle catene di Markov abbiamo una matrice P,
detta matrice di transizione, i cui elementi rappresentano delle
probabilità di transizione tra diversi stati. Più precisamente pij è
la probabilità condizionata che il sistema si trovi “domani” nello
stato j essendo “oggi” nello stato i.
P è una matrice stocastica


la somma di ogni riga di P è uguale ad 1.
Indicato con p0 il vettore iniziale di probabilità degli stati si è
interessati a sapere cosa succede al variare di n al vettore pn
definito come
pn+1 = Ppn = Pnp0, n >= 0
Catene di Markov

Si noti che anche Pn è una matrice di transizione avente righe a
somma 1, un suo generico elemento di indici i e j rappresenta
quindi la probabilità che il sistema si trovi dopo n passi nello
stato j trovandosi nello stato i all’istante iniziale.


elemento ij della matrice Pn
Si è interessati a sapere cosa succede per n crescente. Cosa
possiamo dire sul “comportamento” della matrice Pn e di pn?


un concetto chiave in questo contesto è la distribuzione stazionaria, che
indicheremo d’ora in poi con π.
una distribuzione π si dice stazionaria per una catena con probabilità di
transizione P(x,y) se
p *  p *P
il vettore delle probabilità di essere in un certo stato è indipendente dalla
condizione iniziale.
Esempio

spazio degli stati S={pioggia, sole, nuvole}





 0.5 0.25 0.25 


P   0.5
0
0.5 
 0.25 0.25 0.5 


tra 2 giorni p2= p0P2 =(0.375 0.25 0.375) e tra 7 p2= p0P7 =(0.4 0.2 0.4)
Se oggi piove p0 = (1 0 0)


P(pioggia domani| pioggia oggi) = 0.5
P(sole domani| pioggia oggi) = 0.25
P(nuvole domani| pioggia oggi) = 0.25
Se oggi c’è il sole p0 = (0 1 0)


il tempo segue un processo di Markov: la probabilità del tempo di domani
dipende solo da quello di oggi
tra 2 giorni p2= p0P2 =(0.4375 0.1875 0.375) e tra 7 p2= p0P7 =(0.4 0.2 0.4)
Dopo un certo tempo il tempo atteso è indipendente dal vettore
delle probabilità iniziali

la catena ha raggiunto una distribuzione stazionaria, dove i valori di
probabilità sono indipendenti dai valori iniziali di partenza
Irriducibilità e aperiodicità

Una catena di Markov è irriducibile se
 nij  0 t.c. pij ij  0 i, j
n


in altre parole tutti gli stati comunicano con gli altri, è sempre possibile
muoversi da uno stato all’altro
Una catena di Markov è aperiodica se il numero di step richiesti
per muoversi tra due stati non è un multiplo di qualche intero. La
catena non è intrappolata in qualche ciclo di lunghezza fissata tra
certi stati


il periodo di uno stato x è il massimo comune divisore dell’insieme
dx={n ≥ 1 : Pn(x, x) > 0}.
uno stato x si dice aperiodico se dx = 1.
Ergodicità

Se {Xt} è una catena di Markox aperiodica e irriducibile con
distribuzione stazionaria p

p è l’unica distribuzione invariante

1
f ( X i )   f ( x)p ( x)dx  E{ f ( X )}
lim
N  N


N
P
lim ( x, y)  p ( X )
N 
La media campionaria degli stati della catena è stimatore
consistente del valore atteso della distribuzione limite π,
nonostante gli stati siano fra loro dipendenti.
Catene reversibili

Condizione sufficiente per avere una distribuzione stazionaria è
la condizione di reversibilità
P ( j , k )p *j  p *j P (k , j ) k , j

Una catena di Markov che soddisfa tale condizione si dice
reversibile


L’importanza delle catene reversibili si spiega immediatamente:


tale condizione implica che la catena ammette distribuzione stazionaria
se esiste una distribuzione π che soddisfa la condizione di reversibilità per
una catena di Markov irriducibile e aperiodica la distribuzione π è
distribuzione stazionaria e anche distribuzione limite.
La costruzione di catene di Markov con distribuzione limite si
riduce a trovare probabilità di transizione che soddisfano la
condizione di reversibilità
MCMC

La strategia di campionamento MCMC consiste nella costruzione
di catene di Markov aperiodiche e irriducibili per le quali la
distribuzione stazionaria sia esattamente la distribuzione target
π.

Due sono gli algoritmi utilizzati nel contesto della simulazione
MCMC:


Algoritmo Metropolis-Hasting
Algoritmo Gibbs Sampler
MCMC

Sia S = {s1, s2, . . . , sm} un insieme di stati, dove m è la
cardinalità dell’insieme S, e sia X una variabile aleatoria a valori in
S, con probabilità pj = P(X = sj ). Sia, inoltre, f una funzione
definita su S a valori in R. Si vuole calcolare
m
  E[ f ( X )]   f ( s j ) p j
j 1

Tale quantità potrebbe essere calcolata direttamente utilizzando la formula
sopra riportata. Tuttavia, se la cardinalità di S è molto grande, tale approccio
può risultate eccessivamente oneroso. Alternativamente, la quantità potrebbe
essere approssimata utilizzando il metodo di Monte Carlo, ovvero
campionando la variabile X, e utilizzando lo stimatore di media campionaria.
Questo presuppone, tuttavia, di saper campionare dall’insieme S, cosa non
sempre scontata. Inoltre, in talune applicazioni, la densità di probabilità è nota
soltanto a meno di una costante moltiplicativa il che rende impossibile
l’utilizzo delle tecniche menzionate.
MCMC

L’idea dei metodi Markov Chain Monte Carlo è la seguente se
siamo in grado di costruire una catena di Markov Xn sugli stati S,
che sia ergodica e abbia p come probabilità limite, possiamo
approssimare la quantità  mediante
n 1
n   f ( X i ) n
i 0


Per la proprietà di ergodicità, sappiamo infatti
lim n   n  
che qualunque sia il punto di partenza della catena.
Poiché i primi stati della catena sono fortemente influenzati dallo stato
iniziale, è buona norma iniziare a calcolare la media campionaria, lungo la
traiettoria della catena, soltanto dopo k iterazioni, con k scelto
opportunamente. Definiamo dunque lo stimatore
 k ,n
1 k n
  f (Xi )
n i  k 1
Metropolis-Hastings

Presenteremo ora l’algortimo MCMC universalmente noto con il
nome di Metropolis-Hastings.


Esso risale all’originale idea di Metropolis, alla base di svariati algoritmi di
campionamento (i.e. simulated annealing): l’algoritmo è basato
sull’analogia termodinamica con la posizione di equilibrio di un certo
numero di molecole in una sostanza, la cui distribuzione è data da
un’energia potenziale; dal momento che il campionamento diretto da
questa distribuzione non è possibile, Metropolis propose l’utlizzo di
metodi Monte Carlo.
In seguito, Hastings riprese l’idea originale inserendola nel framework del
campionamento da catene di Markov, e mantenendo l’accettazione o il
rifiuto del valore campionato nel nucleo di transizione della catena.
Metropolis-Hastings



L’idea è la seguente: supponiamo di poter costruire una catena di
Markov {Xn} irriducibile e aperiodica, con un’unica legge limite
p
Se noi simuliamo l’evoluzione della catena la legge di Xn al
tempo n, quando n ->∞, convergerà a p, indipendentemente
dalla legge iniziale scelta.
Consideriamo una matrice stocastica Q  R mxm irriducibile e
aperiodica t.c. qij≠0 se qji≠0 , i,j=1,…,m. Sia Yn la catena di
Markov sugli stati S avente Q come matrice di transizione. La
catena Yn non avrà, in generale, probabilità limite pari a p.
Metropolis-Hastings

Costruiamo la matrice di transizione P definita da
pij   ij qij , i  j
 ij  min(1,

q jip j
qij p i
j i
)
La catena Xn è ergodica, reversibile e ammette p come
distribuzione limite

L’equazione di bilancio è infatti soddisfatta
p i min(1,
q jip j
qij p i
)qij  p j min(1,
analizzando i due casi q jip j  qij p i

pii  1   pij
qij p i
q jip j
)q ji
q jip j  qij p i l’eq è sempre soddisfatta
Dall’equazione di bilancio dettagliato consegue che p è anche
probabilità invariante di P, e quindi probabilità limite, essendo la
catena ergodica:
Metropolis-Hastings - algoritmo

Un modo per generare il nuovo stato Xk+1 della catena a partire
da Xk

campionare la variabile Y avente distribuzione qXkj e calcolare la
probabilità di accettazione dello spostamento da xk a y
 ij  min(1,


q jip j
qij p i
)
accettare y con tale probabilità; come?
estrarre t in [0,1] da una distribuzione di probabilità uniforme


se  ( xk , y)  t , set xk 1  y, else xk 1  xk
se k=K stop else k=k+1 and go to step 2
Esercizio



Si consideri un sistema di particelle che possa assumere solo m
configurazioni possibili S = {1, 2, . . . ,m}. La probabilità che il
sistema si trovi nella configurazione j-esima è pj = Ce−Ej/T
(distribuzione di Boltzmann), essendo Ej = j2 l l’energia del
sistema, T la temperatura e C la costante di normalizzazione.
Scrivere una funzione Matlab che implementi l’algoritmo di
Hastings-Metropolis, dati i valori m e T, il numero n di passi da
simulare e lo stato iniziale X0 della catena e restituisca il vettore
X degli stati visitati. Si prenda la matrice di transizione qij = 1/m
(ovvero partendo dallo stato i, lo stato candidato è uno
qualunque degli altri stati con uguale probabilità).
Si prenda m = 100, T = 100 e si parta da uno stato a caso. Si
n
valuti lo stimatore. k , n   E X4 / n
confrontandolo col
i  k 1
valore esatto
i
Esercizio

S = {1, 2, . . . ,m},
pj = Ce−Ej/T , Ej = j2
function X = my_mh(m,T,n,X0)


n numero di passi da simulare, stato iniziale X0 della catena , X vettore
degli stati visitati, matrice di transizione qij = 1/m (partendo dallo stato i,
lo stato candidato è uno qualunque degli altri stati con uguale probabilità)
Algoritmo



campionare la variabile Y avente distribuzione qXkj e calcolare la
probabilità di accettazione dello spostamento da xk a y. P(Xk=si)=pi
q jip j
 ij  min(1,
)
qij p i
accettare y con tale probabilità; come?
estrarre u in [0,1] da una distribuzione di probabilità uniforme
 ij  u, set

xk 1  y, else xk 1  xk
se k=K stop else k=k+1 and go to step 2
Metropolis-Hastings




Costruiamo un opportuno nucleo di transizione P(x, y) tale che p
sia la distribuzione stazionaria
Un modo immediato è scegliere una q tale che:
p(x)P(x, y) = p(y)P(y, x), ∀ (x, y)
In tal caso la catena è reversibile, condizione sufficiente perché p
sia distribuzione stazionaria della catena risultante.
Il nucleo P(x, y) è scelto tale che
P(x, y) = q(x, y)α(x, y), se x !=y,


dove q(x, y) è un nucleo di transizione arbitrario, mentre α(x, y) è una
probabilità di accettazione.
la probabilità di accettazione α(·, ·) è scelta in modo che la catena
risultante sia reversibile, ovvero
q ( y, x)p ( y )
 ( x, y )  min(1,
)
q ( x, y )p ( x)
Metropolis-Hastings



Per dimostrare che tale algoritmo genera una catena di Markov la
cui densità di equilibrio è p(x) è sufficiente mostrare che il
transition kernel P soddisfa l’equazione di bilancio
Noi campioniamo da q(x,y) e accettiamo di muoverci con
probabilità α(x, y)
Il transition kernel è
P(x,y)=q(x,y) α(x, y)


si dimostra che con questo transition è soddisfatta l’equazione di bilancio
la dimostrazione che la condizione di reversibilità è soddisfatta dalla scelta di P e
dunque definisce una catena reversibile con distribuzione di equilibrio π, segue
immediatamente dalla definizione della probabilità di accettazione.
Metropolis-Hastings

La scelta del nucleo q(·, ·) è arbitraria, e fornisce uno strumento
molto flessibile per la costruzione dell’algoritmo;
5
2
1
3
4
1
1
2
3
4
4
5
…
Metropolis-Hastings

Il nucleo di transizione q definisce solo una possibile mossa della
catena, che deve essere confermata in base ad α; per tale ragione
viene solitamente chiamato proposal.


La catena potrebbe rimanere nello stesso stato per molte
iterazioni: la potenza del metodo si esplica quando si riesce ad
avere un tasso di accettazione non troppo basso.


q deve dosare opportunamente i movimenti proposti in modo che non
risulti troppo basso, ma lo spazio degli stati venga visitato sufficientemente
in fretta.
Monitoraggio: percentuale di iterazioni in cui la proposta è accettata.
Deve essere facile campionare dalla proposal, in quanto il
metodo sostituisce il campionamento da p (difficile) con molti
campionamenti da q (facili);
Metropolis-Hastings - algoritmo


inizializzare il contatore delle iterazioni k=0 ed il valore iniziale
x0 per lo stato della catena
estrarre y da una proposal distribution q(xk,y) e calcolare la
probabilità di accettazione dello spostamento da xk a y
 ( xk , y )  min(1,

q ( y, xk )p ( y )
)
q ( xk , y )p ( xk )
accettare y con tale probabilità; come?
estrarre t in [0,1] da una distribuzione di probabilità uniforme


se  ( xk , y)  t , set xk 1  y, else xk 1  xk
se k=K stop else k=k+1 and go to step 2
Metropolis

L'algoritmo base genera una sequenza stocastica di stati, a partire
da x0, e la seguente regola dinamica per passare dallo stato x a y:

propone un candidato y a partire da q(y,x) che potrebbe dipendere da x. Se
il kernel è simmetrico q(x,y)=q(y,x) per ogni x,y allora
 ( x, y )  min(1,

p ( y)
)
p ( x)
La generazione del candidato si può implementare come y = x +
e.


Possibili scelte per la distribuzione di e sono distribuzioni di Cauchy,
Gaussiane, o T con pochi gradi di liberta.
se y usa meno energia dello stato corrente x si accetta con probabilità 1.
Altrimenti si accetta con una probabilità esponenzialmente decrescente
nella dierenza di energia; formalmente


inizializzare il contatore delle iterazioni k=0 ed il valore iniziale
x0 per lo stato della catena
estrarre y da una proposal distribution q(xk,y) e calcolare la
probabilità di accettazione dello spostamento da xk a y
 ( xk , y )  min(1,

q( y, xk )p ( y )
)
q( xk , y )p ( xk )
accettare y con tale probabilità; come?
estrarre t in [0,1] da una distribuzione di probabilità uniforme


se  ( xk , y)  t , set xk 1  y, else xk 1  xk
se k=K stop else k=k+1 and go to step 2
1
p ( x) 
1 x2
q ( x, y )  exp(
1
2 2
x y )
2
Esercizio

Si assuma di voler campionare da una densità di Cauchy non
normalizzata
1
y
1 x2

Usiamo il random walk come proposal
q( x, y )  exp(
1
2 2
x y )
2
Riassumendo



Le tecniche Monte Carlo per l’approssimazione di integrali
possono essere contestualizzate in ambito statistico (inferenziale)
e sfruttate per il calcolo di tutte quelle quantità di interesse
statistico esprimibili sottoforma di integrali.
Le tecniche Markov Chain Monte Carlo, che permettono di
ottenere un campione da una qualsiasi distribuzione target di
interesse, consentono tra le altre cose di utilizzarlo per calcolare
stime e fare inferenza. In questo senso possono considerarsi una
generalizzazione delle tecniche Monte Carlo.
I metodi di inferenza statistica Bayesiana, basati sulla simulazione
stocastica, utilizzano campioni dalla distribuzione a posteriori
(target) per riassumere l’informazione in essa contenuta.