Laboratorio Processi Stocastici

1
Laboratorio Processi Stocastici

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

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

3
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)
• La velocità di convergenza dipende da N e non
  dallo spazio nel quale vivono i campioni
• principale vantaggio rispetto alle tecniche di
  quadratura

4
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(xy), 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
  gt per ottenere un insieme di campioni con cui
  approssimare lintegrale I(f) si ricorre a
  tecniche più sofisticate come limportance
  sampling o i metoti MCMC

5
Importance sampling

• Lidea alla base dellIS è 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 random sampling
  di f(x)w(x)

6
Importance sampling

• Se sono in un contesto di inferenza Bayesiana
  vorrei poter campionare p(xy)
• campionare tale pdf usando lIS 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 nellintegrazione
  Monte Carlo.
• Lalgoritmo Metropolis è considerato tra uno dei
  dieci metodi che hanno avuto maggiore influenza
  sulle scienze e lingegneria nel XX secolo

7
MCMC
8
Da MC a MCMC

• Lobiettivo dei metodi MC, e il loro principale
  utilizzo nelle applicazioni, è lintegrazione
• 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.

9
MCMC

• Lidea di base dei metodi MCMC è di generare un
  campione dalla distribuzione dinteresse 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 .
• Lapplicazione più popolare dei metodi MCMC è
  nellambito dellinferenza Bayesiana, dove la
  distribuzione target p è la distribuzione a
  posteriori, generalmente indicata con p, ed X
  sono i parametri di interesse, generalmente
  indicati con ?.

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

11
Catene di Markov

• Sia Xt il valore di una v.a. a t e sia Ss1, ,
  sm linsieme degli stati. Il processo stocastico
  Xt è un processo di Markov se

12
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
• pn1 Ppn Pnp0, n gt 0

13
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 allistante 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 dora
  in poi con p.
• una distribuzione p 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.

14
Esempio

• spazio degli stati Spioggia, 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

15
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 allaltro
• 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 dellinsieme dxn 1 Pn(x,
  x) gt 0.
• uno stato x si dice aperiodico se dx 1.

16
Ergodicità

• Se Xt è una catena di Markox aperiodica e
  irriducibile con distribuzione stazionaria p
• p è lunica distribuzione invariante
• La media campionaria degli stati della catena è
  stimatore consistente del valore atteso della
  distribuzione limite p, nonostante gli stati
  siano fra loro dipendenti.

17
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
• Limportanza delle catene reversibili si spiega
  immediatamente
• se esiste una distribuzione p che soddisfa la
  condizione di reversibilità per una catena di
  Markov irriducibile e aperiodica la distribuzione
  p è 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à

18
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 p.
• Due sono gli algoritmi utilizzati nel contesto
  della simulazione MCMC
• Algoritmo Metropolis-Hasting
• Algoritmo Gibbs Sampler

19
MCMC

• Sia S s1, s2, . . . , sm un insieme di stati,
  dove m è la cardinalità dellinsieme 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
  dallinsieme 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
  lutilizzo delle tecniche menzionate.

20
MCMC

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

21
Metropolis-Hastings

• Presenteremo ora lalgortimo MCMC universalmente
  noto con il nome di Metropolis-Hastings.
• Esso risale alloriginale idea di Metropolis,
  alla base di svariati algoritmi di campionamento
  (i.e. simulated annealing) lalgoritmo è basato
  sullanalogia termodinamica con la posizione di
  equilibrio di un certo numero di molecole in una
  sostanza, la cui distribuzione è data da
  unenergia potenziale dal momento che il
  campionamento diretto da questa distribuzione non
  è possibile, Metropolis propose lutlizzo di
  metodi Monte Carlo.
• In seguito, Hastings riprese lidea originale
  inserendola nel framework del campionamento da
  catene di Markov, e mantenendo laccettazione o
  il rifiuto del valore campionato nel nucleo di
  transizione della catena.

22
Metropolis-Hastings

• Lidea è la seguente supponiamo di poter
  costruire una catena di Markov Xn irriducibile
  e aperiodica, con ununica legge limite p
• Se noi simuliamo levoluzione della catena la
  legge di Xn al tempo n, quando n -gt8, convergerà
  a p, indipendentemente dalla legge iniziale
  scelta.
• Consideriamo una matrice stocastica
  irriducibile e aperiodica t.c. qij?0 se qji?0 ,
  i,j1,,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.

23
Metropolis-Hastings

• Costruiamo la matrice di transizione P definita
  da
• La catena Xn è ergodica, reversibile e ammette p
  come distribuzione limite
• Lequazione di bilancio è infatti soddisfatta
• analizzando i due casi
  leq è sempre soddisfatta
• Dallequazione di bilancio dettagliato consegue
  che p è anche probabilità invariante di P, e
  quindi probabilità limite, essendo la catena
  ergodica

24
Metropolis-Hastings - algoritmo

• Un modo per generare il nuovo stato Xk1 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
• accettare y con tale probabilità come?
• estrarre t in 0,1 da una distribuzione di
  probabilità uniforme
• se
• se kK stop else kk1 and go to step 2

25
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
  lenergia del sistema, T la temperatura e C la
  costante di normalizzazione.
• Scrivere una funzione Matlab che implementi
  lalgoritmo 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

26
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(Xksi)pi
• accettare y con tale probabilità come?
• estrarre u in 0,1 da una distribuzione di
  probabilità uniforme
• se kK stop else kk1 and go to step 2

27
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)a(x, y), se x !y,
• dove q(x, y) è un nucleo di transizione
  arbitrario, mentre a(x, y) è una probabilità di
  accettazione.
• la probabilità di accettazione a(, ) è scelta
  in modo che la catena risultante sia reversibile,
  ovvero

28
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 lequazione di bilancio
• Noi campioniamo da q(x,y) e accettiamo di
  muoverci con probabilità a(x, y)
• Il transition kernel è
• P(x,y)q(x,y) a(x, y)
• si dimostra che con questo transition è
  soddisfatta lequazione di bilancio
• la dimostrazione che la condizione di
  reversibilità è soddisfatta dalla scelta di P e
  dunque definisce una catena reversibile con
  distribuzione di equilibrio p, segue
  immediatamente dalla definizione della
  probabilità di accettazione.

29
Metropolis-Hastings

• La scelta del nucleo q(, ) è arbitraria, e
  fornisce uno strumento molto flessibile per la
  costruzione dellalgoritmo

5
2
1
3
4

1
2
3
4
5
1
4
30
Metropolis-Hastings

• Il nucleo di transizione q definisce solo una
  possibile mossa della catena, che deve essere
  confermata in base ad a 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)

31
Metropolis-Hastings - algoritmo

• inizializzare il contatore delle iterazioni k0
  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
• accettare y con tale probabilità come?
• estrarre t in 0,1 da una distribuzione di
  probabilità uniforme
• se
• se kK stop else kk1 and go to step 2

32
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

33

• inizializzare il contatore delle iterazioni k0
  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
• accettare y con tale probabilità come?
• estrarre t in 0,1 da una distribuzione di
  probabilità uniforme
• se
• se kK stop else kk1 and go to step 2

34
Esercizio

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

35
Riassumendo

• Le tecniche Monte Carlo per lapprossimazione 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 linformazione in essa
  contenuta.