ML - 15 - Montecarlo Methods

Lecture Info

Date: [2021-04-04 dom 23:40]
Lecturer: Giorgio Gambosi
Slides: (ml_07_sampling-slides.pdf)

In questa lezione abbiamo introdotto i metodi di montecarlo, con una particolare enfasi sull'approccio generale Markov Chain Montecarlo (MCMC) che ci permette di effettuare del sampling da una generica distribuzione $p(x)$ a patto di poter calcolare un valore $\pi(x) = K \cdot p(x)$ che sia proporzionale a meno di una qualche costante ignota $K$. Dopo aver fatto vedere qualche variante MCMC abbiamo terminato la lezione discutendo come applicare queste tecniche nel contesto dell'apprendimento bayesiano.

1 The Basic Problem

I metodi montecarlo sono dei metodi numerici che offrono svariate applicazioni in vari contesti. Inizialmente sono stati introdotti per approssimare numericamente degli integrali la cui risoluzione attraverso metodi analitici risulta essere troppo difficile, sia per la particolare funzione da integrare e sia per l'eventuale elevata dimensionalità della stessa.

\[\int_a^b g(x) \,\, dx\]

L'idea di base è quella di vedere l'integrale come calcolo di un valore atteso. In particolare assumiamo di avere una funzione $f(x)$ e una densità di probabilità $p(x)$ entrambi definiti nell'intervallo $[a, b]$ e tali che $g(x) = f(x) \cdot p(x)$. Allora possiamo scrivere

\[\int_a^b g(x) \,\, dx = \int_a^b f(x) \cdot g(x) \,\, dx = \mathbb{E}_{p(x)}[f(x)]\]

l'idea è quindi di approsimare questo valore atteso utilizzando $n$ valori $f(x_1), \ldots, f(x_n)$ dove i punti $x_1, \ldots, x_n$ sono ottenuti campionando $p(x)$

\[\mathbb{E}[f(x)] \approx \frac{1}{n} \sum\limits_{i = 1}^n f(x_i)\]

L'approccio utilizzato è quindi il seguente:

Si effettua il campionamento di una sequenza di valori $x^{(1)}, x^{(2)}, \ldots, x^{(n)}$ dalla distribuzione $p(x)$, ovvero tali che $P(X = x^{(i)}) = p(x^{(i)})$.
Si applica la funzione $f(x)$ ai valori campionati.
Si calcola la media dell'insieme di valori ottenuti.

Questo ovviamente comporta di riuscire a ottenere dei campioni dalla distribuzione $p(x)$, che essendo una distribuzione generica potrebbe essere difficile da descrivere in modo analitico.

Per risolvere questo problema si effettua la seguente ipotesi: si assume di avere un generatore (pseudo) casuale $\mathcal{R}$ che ritorna una sequenza di valori che sono (approssimatamente) distribuiti in modo uniforme nell'intervallo $[0, 1]$. Data la distribuzione $p(x)$, la sequenza di valori ritornata da $\mathcal{R}$ può poi essere utilizzata per derivare una diversa sequenza di valori distribuiti secondo $p(x)$.

2 Sampling General Distributions

Così facendo abbiamo quindi ridotto il problema di stimare numericamente il valore di un integrale definito siamo adesso giunti al problema di come effettuare del sampling da una generica distribuzione, specialmente nelle situazioni in cui non abbiamo a disposizione nessun tipo di descrizione analitica di tale distribuzione.

A parte il caso semplice in cui conosciamo l'espressione analitica di $p(x)$ infatti, nella maggior parte dei casi non è possibile derivare immediatamente valori distribuiti secondo $p(x)$. A tale fine nel corso degli anni sono stati introdotto vari metodi per effettuare questo tipo di sampling:

rejection sampling.
importance sampling.
adaptive rejection sampling.
sampling-importance-sampling.

2.1 Easy Case

Il caso più semplice è quando conosciamo la descrizione analitica della distribuzione $p(x)$ e siamo interessati a trovare una funzione $f(z)$ tale che se $z \sim U(0, 1)$, allora $f(z) \sim p(z)$.

Utilizzando il linguaggio delle funzioni di probabilità cumulative la nostra richiesta è equivalente a dire che per ogni $z \in [0,1]$, la probabilità cumulativa di $z$ rispetto alla distribuzione uniforme deve essere uguale alla probabilità cumulativa di $f(z)$ rispetto a $p(z)$. In formula

\[z = \int\limits_0^z 1 \,\, d\zeta = \int\limits_{-\infty}^{f(z)} p(\zeta) \,\, d\zeta = P(f(z))\]

da questo deriva che

\[z = P(f(z)) \iff f(z) = P^{-1}(z)\]

dove $P(\cdot)$ è la distribuzione cumulativa di $p(x)$.

2.1.1 Example 1: Exponential

Come primo esempio assumiamo che il nostro PRNG $\mathcal{R}$ ci permette di produrre valori che seguono una distribuzione esponenziale, ovvero tali che

\[P(x|\lambda) = \lambda e^{-\lambda x} \,\,\,,\,\,\, 0 \leq x < \infty\]

dove la distribuzione cumulativa di una esponenziale è descritta come segue

\[P(x) = \int\limits_0^x \lambda e^{- \lambda \xi} \,\, d \xi = 1 - e^{-\lambda x}\]

In questo caso quindi impostando $z = P(f(z)) = 1 - e^{- \lambda f(z)}$ otteniamo

\[\begin{split} & \,\,\, e^{- \lambda f(z)} = 1 - z \\ \iff \,\,\, & - \lambda f(z) = \ln{(1 - z)} \\ \iff \,\,\, & f(z) = -\frac{1}{\lambda} \ln{(1 - z)} \\ \end{split}\]

Questo significa che applicando la funzione $f(z)$ agli elementi generati dal nostro PRNG $\mathcal{R}$ otteniamo una nuova sequenza distribuita secondo la distribuzione $\lambda e^{-\lambda x}$.

2.2 Rejection Sampling

Nel rejection sampling assumiamo che $p(x)$ sia difficile da campionare in quanto non abbiamo nessuna descrizione analitica di tale distribuzione, ma che esista una distribuzione $q(x)$ più facile da campionare e tale che esiste un $K$ tale che $K q(x) \geq p(x)$ per ogni $x$.

Il metodo è così descritto:

Si campiona $\overline{x} \sim q(x)$.
Si campiona $u^* \sim U(0, Kq(\overline{x}))$.
Se $u^* \leq p(\overline{x})$ si accetta e si ritorna l'elemento campionato, altrimenti lo si scarta.

2.3 Importance Sampling

In questo caso assumiamo di voler approssimare il valore atteso di $f(x)$ rispetto a $p(x)$

\[\mathbb{E}_{p(x)}[f(x)] = \int_x f(x) \cdot p(x) \,\, dx \]

Si suppone quindi che $p(x)$ sia difficile da campionare ma sia facile da calcolare per ogni punto $x$, e sia $q(x)$ una distribuzione che può essere facilmente campionata. Allora,

\[\int_x f(x) p(x) \,\,\, dx = \int_x \frac{f(x)}{q(x)} p(x) q(x) \,\,\, dx\]

A questo punto cambiamo prospettiva e vediamo $q(x)$ come la funzione da campionare e $f(x)/q(x) \cdot p(x)$ come la funzione da applicare. Così facendo l'approssimazione è calcolata estraendo una sequenza $x^{(1)}, \ldots, x^{(n)}$ di $n$ punti che seguono una distribuzione $q(x)$, e calcolando per ogni punto il valore $f(x)/q(x) \cdot p(x)$, ottenendo

\[\int_x \frac{f(x)}{q(x)} p(x) q(x) \,\,\, dx \approx \frac{1}{n} \sum\limits_{i = 1}^n \frac{p(x^{(i)})}{q(x^{(i)})} \cdot f(x^{(i)})\]

3 Markov Chain Montecarlo

Uno strumento estremamente utile per affrontare il problema di generare degli elementi $x_1, x_2, \ldots, x_n$ che seguono una determinata distribuzione $p(x)$ ci viene offerto dalle Markov chain (catene di Markov). Questo metodo prende il nome di MCMC o Markov Chain MonteCarlo in quanto combina l'estrazione casuale di valori da qualche distribuzione (montecarlo) utilizzando una catena di markov (markov chain).

In un approccio MCMC l'ipotesi che facciamo è che sia possibile valutare la distribuzione $p(x)$. In realtà i metodi MCMC richiedono un qualcosa di più debole: ci basta infatti che sia possibile valutare il valore di $p(x)$ a meno di una costante, ovvero ci basta essere in grado di calcolare

\[\pi(x) = K \cdot p(x)\]

dove la costante $K$ può rimanere ignota.

3.1 Markov Chains

Ricordiamo brevemente che, data una sequenza (possibilmente infinita) di variabili aleatorie $\mathbf{X} = (X_0, X_1, \ldots)$ e uno spazio degli stati $\chi$ di tutti i possibili valori per ogni v.a. $X_i \in \mathbf{X}$, una catena di Markov su $\mathbf{X}$ è un processo stocastico che definisce, per ogni coppia ordinata $(x_i, x_j) \in \chi^2$, una probabilità $p_{i, j}$ di transizione da $x_i$ a $x_j$ tale che per ogni $t > 0$

\[P(X_t = x_i | X_{t-1} = x_j) = p_{i,j}\]

Dato uno stato iniziale, ovvero un valore per $X_0$, la distribuzione $P(X_t = x_j | X_0 = x_k)$ per ogni variabile aleatoria $X_t$ nell'insieme degli stati può essere facilmente calcolata tramite moltiplicazioni matriciali.

È possibile dimostrare che sotto opportune condizioni una catena di Markov è detta ergodica, ovvero che la probabilità $P(X_t = x_j | X_0 = x_k)$ per $n \to \infty$ è

Indipendente dallo stato iniziale

\[P(X_t = x_j | X_0 = x_k) = P(X_t = x_j)\]
Stazionaria

\[P(X_t = x_j) = P(X_{t+1} = x_j)\]

3.2 MCMC Idea

Quando abbiamo a che fare con una catena di Markov ergodica è possibile stimare la distribuzione stazionaria andando a campionare la frequenza con cui la catena entra nei vari stati. In particolare si potrebbe aspettare un po' di tempo, ad esempio $n = 10000$ transizioni, e poi iniziare a contare quante volte la catena entra in ogni stato per andare a campionare la sua distribuzione stazionaria.

Quindi, se riuscissimo a trovare una catena di markov la cui distribuzione stazionaria è proprio $p(x)$, allora procedendo come descritto sopra saremmo anche in grado di campionare $p(x)$.

Per effetuare questo ovviamente necessitiamo che l'insieme di stati della catena di markov deve corrispondere al dominio della distribuzione di interesse. In particolare se $p(x)$ è definita sui reali $\mathbb{R}$ questo significa che dovremmo lavorare con catene di markov con un numero reale di stati.

Data una distribuzione difficile da campionare $p(x)$ l'idea dietro a MCMC è quindi quella di costruire una catena di Markov ergodica tale che

Le probabilità di transizione $P(x_i | x_{i-1})$ sono facili da campionare.
La distribuzione stazionaria è $p(x)$.

3.3 How to Use it

Data una catena di Markov per applicare un metodo di tipo MCMC si procede come segue:

Come prima cosa si effettua una sequenza di un certo numero di transizioni random partendo da un qualsiasi stato iniziale. Questa sequenza di transizioni è chiamata burn-in time.
Dopo questo numero di transizioni iniziali ad ogni step il valore $\overline{x}$ raggiunto dalla catena di Markov viene testato rispetto a un criterio predefinito: se il test è positivo, il valore è ritornato.

I valori ritornati sono (approssimatamente) distribuiti come $p(x)$, e quindi la sequenza di valori ritornati può essere utilizzata come sequenza campionata da $p(x)$.

Vari MCMC metodi sono stati definiti nel corso del tempo, ciascuno dei quali definisce una propria catena di Markov con una particolare struttura con dei propri criteri di accettazione.

4 Metropolis Algorithm

Sia $x^{(i-1)}$ lo stato della catena di Markov subito dopo il burning time, e sia $\overline{x}$ il valore prodotto da una transizione random a partire da $x^{(i-1)}$ ottenuta effettuando il sampling di $q(x | x^{(i-1)})$, che è proprio la distribuzione di probabilità degli stati di essere raggiunti da una transizione a partire dallo stato $x^{(i-1)}$.

A questo punto accettiamo il valore $\overline{x}$ e lo ritorniamo come campione di $p(x)$ con la seguente probabilità

\[A(\overline{x}, x^{(i-1)}) = \min{\Bigg(1, \,\, \frac{p(\overline{x})}{p(x^{(i-1)})}\Bigg)}\]

Nel caso in cui $\overline{x}$ è accetato, allora $x^{(i)} = \overline{x}$ diventa lo stato corrente della catena di Markov, altrimenti lo stato corrente non viene modificato, ovvero $x^{(i)} = x^{(i-1)}$.

La struttura dell'algoritmo di Metropolis è quindi la seguente:

Estraggo un valore $\overline{x}$ da $q(x | x^{(i-1)})$.
Con probabilità $A(\overline{x}, x^{(i-1)})$
- definisco $x^{(i)} = \overline{x}$ e ritorno $\overline{x}$.
- altrimenti definisco $x^{(i)} = x^{(i-1)}$.

È possibile dimostrare che per ogni coppia $x, x^{'}$ la probabilità reale di transizioni da $x$ a $x^{'}$ è data dalla transition kernel

\[T(x|x^{'}) = q(x|x^{'}) \cdot A(x|x^{'})\]

dove $q(x|x^{'})$ è la probabilità di transizione e $A(x|x^{'})$ è la probabilità di accettazione. Così facendo quindi stiamo costruendo una nuova catena di Markov la cui distribuzione stazionaria è proprio $p(x)$.

4.1 Existence of Stationary Distribution

Data una distribuzione $p(x)$, una catena di Markov ha la proprietà detta detailed balance rispetto a $p(x)$ se, per ogni coppia $x, x^{'}$, si ha che

\[p(x) \cdot q(x^{'} | x) = p(x^{'}) q(x | x^{'})\]

ovvero che la probabilità che ad un certo punto lo stato corrente è $x$ e il successivo è $x^{'}$ è uguale alla probabilità che lo stato corrente è $x^{'}$ e che il prossimo è $x$.

È possibile dimostrare che se questa proprietà vale, allora $p(x)$ è proprio la distribuzione stazionaria della catena di Markov, in quanto se $p^{*}(x)$ è la distribuzione stazionaria, allora per definizione

\[p^{*}(x) = \sum\limits_{x^{'}} q(x|x^{'}) p^{*}(x^{'})\]

e per $p(x)$ otteniamo

\[p(x) = \sum\limits_{x^{'}} q(x|x^{'})p(x^{'}) = \sum\limits_{x'}q(x^{'}|x)p(x) = p(x) \sum\limits_{x^{'}}q(x^{'}|x) = p(x)\]

4.2 Uniqueness of Stationary Distribution

Anche nel caso in cui $p(x)$ è la distribuzione stazionaria della nostra catena, per essere sicuri che la catena tende a $p(x)$ per ogni stato iniziale dobbiamo verificare che la catena si ergodica.

Una condizione sufficiente a tale fine è controllare che per ogni coppia $x, x^{'}$ la probabilità di transizione è positiva, ovvero $q(x | x^{'}) > 0$.

4.3 Why it Works

Assumendo di avere a disposizione una distribuzione delle probabilità di transizione $q(\cdot)$ che sia

symmetric: $q(x|x^{'}) = q(x^{'}|x)$
positive: $q(x|x^{'}) > 0$

allora per far vedere che l'algoritmo di Metropolis funzione ci basta semlicemente far vedere che la balanced property vale rispetto a $p(x)$ per la catena di Markov descritta dalle seguenti probabilità di transizione

\[T(x|x^{'}) = q(x|x^{'}) \cdot A(x, x^{'})\]

Questo può essere dimostrato nel seguente modo

\[\begin{split} p(x) \cdot T(x^{'}|x) = \ldots = p(x^{'}) \cdot T(x|x^{'}) \end{split}\]

\[\tag*{$\blacksquare$}\]

Dato che tutte le probabilità di transizione sono positive, otteniamo che la nuova catena di markov è ergodica con distribuzione stazinaria $p(x)$.

5 Metropolis-Hastings Algorithm

Nell'algoritmo di metropolis si assumeva che la probabilità di transizione $q(\cdot)$ era simmetrica. Nel caso in cui questo non vale possiamo applicare l'algoritmo di Metropolis-Hastings, in cui la transition kernel è data da

\[T^{'}(x | x^{'}) = q(x | x^{'}) \cdot A^{'}(x, x^{'})\]

dove

\[A^{'}(x, x^{'}) = \min{\Bigg( 1, \,\, \frac{p(x) q(x^{'} | x)}{p(x^{'}) q(x | x^{'})}\Bigg)}\]

Dato che la proprietà di detailed balance continua a valere, otteniamo comunque un algoritmo funzionante per campionare i valori di $p(x)$.

6 Gibbs Sampling

Utilizzando nel caso in cui devo generare dei vettori, questo metodo di sampling...

7 MCMC and Bayesian Models

Andiamo adesso a vedere come poter utilizzare queste tecniche MCMC nel contesto di un apprendimento bayesiano. Ricordiamo a tale fine che in un approccio bayesiano la probabilità a posteriori ha la seguente forma

\[P(\theta | \mathbf{X}) = \frac{P(\mathbf{X} | \theta) P(\theta)}{P(\mathbf{X})} = \frac{P(\mathbf{X} | \theta) P(\theta)}{Z}\]

dove tipicamente $Z$ è difficile da calcolare.

Da quanto discusso in questa lezione, abbiamo visto che i metodi MCMC ci permettono di campionare una distribuzione $p(x)$ assumendo di poter calcolare una funzione proporzionale a $\pi(x) = K p(x)$ per qualche costante $K$.

Questo significa che se siamo in grado di calcolare la prior $P(\theta)$ e la likelihood $P(\mathbf{X}|\theta) = \prod_i P(x_i | \theta)$ per ogni valore di $\theta$, allora possiamo applicare uno dei metodi MCMC visti per effettuare del sampling di $P(\theta | \mathbf{X})$.

7.1 Sampling the Evidence

Notiamo che volendo è possibile applicare quanto detto inizialmente per valutare esplicitamente l'evidenza

\[P(\mathbf{X}) = \int P(\mathbf{X} | \theta) \cdot P(\theta) \,\,\, d \theta\]

come media di $m$ valori $P(\mathbf{X} | \theta_i)$ calcolati campionando $\theta_1, \ldots, \theta_m$ da $P(\theta)$.

\[P(\mathbf{X}) = \int P(\mathbf{X} | \theta) \cdot P(\theta) \,\,\, d \theta \approx \frac{1}{n} \sum\limits_{i = 1}^m P(\mathbf{X} | \theta_i)\]

7.2 Sampling the Predictive Distribution

Le stesse considerazioni di prima possono essere fatte per la distribuzione predittiva

\[P(\mathbf{x} | \mathbf{X}) = \int_{\theta} P(\mathbf{x} | \theta) \cdot P(\theta | \mathbf{X}) \,\,\, d\theta\]

ovvero possiamo stimare la distribuzione predittiva come media di $m$ valori $P(\mathbf{x} | \theta_i)$ ottenuti effettuando un campionamento di $\theta_1, \ldots, \theta_m$ dalla distribuzione a posteriori $P(\theta | \mathbf{X})$.