Laboratorio Processi Stocastici

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

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 -> uso una pdf dalla quale so campionare • MCMC

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) • 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 • 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 è • 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’importancesampling o i metotiMCMC

Importancesampling • L’idea alla base dell’IS è usare una pdf q(x) dalla quale si sa campionare • 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 • in pratica stiamo effettuando un randomsampling di f(x)w(x)

Importancesampling • Se sono in un contesto di inferenza Bayesiana vorrei poter campionare p(x|y) • campionare tale pdf usando l’IS non è sempre possibile in quanto spesso non si può calcolare la costante di normalizzazione • Una possibile soluzione • MarkoxChain 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 • SiaXtil 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

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 Pne 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 il vettore delle probabilità di essere in un certo stato è indipendente dalla condizione iniziale.

Esempio • spazio degli stati S={pioggia, sole, nuvole} • il tempo segue un processo di Markov: la probabilità del tempo di domani dipende solo da quello di oggi • 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) • 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) • 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 • 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 • 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à • Una catena di Markov che soddisfa tale condizione si dice reversibile • tale condizione implica che la catena ammette distribuzione stazionaria • L’importanza delle catene reversibili si spiega immediatamente: • 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 GibbsSampler

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 • 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à q mediante • Per la proprietà di ergodicità, sappiamo infatti 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

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. simulatedannealing): 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 irriducibile e aperiodicat.c.qij≠0 se qji≠0 , i,j=1,…,m. Sia Ynla 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 • La catena Xn è ergodica, reversibile e ammette p come distribuzione limite • L’equazione di bilancio è infatti soddisfatta analizzando i due casi 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 xka y • accettare y con tale probabilità; come? • estrarre t in [0,1] da una distribuzione di probabilità uniforme • se • 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 valuti lo stimatore. confrontandolo col valore esatto

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 xka y. P(Xk=si)=pi • accettare y con tale probabilità; come? • estrarre u in [0,1] da una distribuzione di probabilità uniforme • 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

Metropolis-Hastings • Per dimostrare che tale algoritmo genera una catena di Markov la cui densità di equilibrio è p(x) è sufficiente mostrare che il transitionkernel P soddisfa l’equazione di bilancio • Noi campioniamo da q(x,y) e accettiamo di muoverci con probabilità α(x, y) • Il transitionkernel è 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. • q deve dosare opportunamente i movimenti proposti in modo che non risulti troppo basso, ma lo spazio degli stati venga visitato sufficientemente in fretta. • 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. • 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 proposaldistribution q(xk,y) e calcolare la probabilità di accettazione dello spostamento da xka y • accettare y con tale probabilità; come? estrarre t in [0,1] da una distribuzione di probabilità uniforme • se • 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 • 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 proposaldistribution q(xk,y) e calcolare la probabilità di accettazione dello spostamento da xka y • accettare y con tale probabilità; come? estrarre t in [0,1] da una distribuzione di probabilità uniforme • se • se k=K stop else k=k+1 and go to step 2

Esercizio • Si assuma di voler campionare da una densità di Cauchy non normalizzata • Usiamo il randomwalk come proposal

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.

Laboratorio Processi Stocastici