Sintesi della Posterior

Il problema che ci si trova a dover affrontare in questo articolo riguarda la sintesi della posterior, esponiamo il problema nel modo più generale possibile:

\[E^{\pi(.|x)}[g(\theta)]= \int_{\Theta}{g(\theta) \space \pi(\theta | x)} \space \space d\theta\]

I principali metodi che si possono utilizzare per risolvere il problema appena esposto si possono dividere in tre categorie:

Metodi di integrazione numerica

Metodi analitici

Metodi simulativi

Metodi di integrazione numerica degli integrali

L’integrazione numerica consiste in una serie di metodi che stimano il valore di un integrale definito, senza dover calcolare la primitiva della funzione integranda.

La necessità di utilizzare l’integrazione numerica, è necessaria in quanto non tutte le funzioni ammettono una primitiva in forma esplicita (per esempio la curva gaussiana) oppure la primitiva della funzione può essere molto complicata da valutare o si è davanti ad una funzione disponibile solo in alcuni punti.

I metodi di integrazione numerica possono essere distinti in due macrocategorie:

formule di Newton-Cotes (per esempio metodo dei trapezi);
formule di Gauss (per esempio metodo monte carlo).

Questi metodi per l’approssimazione della posterior funzionano bene solo quando il parametro è molto piccolo, in questo ambito non si utilizzano molto, anche se nel caso di posterior coniugate i metodi monte carlo risultano unottima scelta.

Inoltre spesso non si conosce $m(x)= \int_{\Theta}{f(\underline{x} | \theta) \pi(\theta)} \space \space d\theta$ e questo obbliga ad utilizzare delle approssimazioni.

Metodi analitici

Approssimazione Normale

L’approssimazione Normale è un metodo che utilizza le proprietà asintotiche della distribuzione Normale per riuscire a risolvere la problematica della sintesi della posterior quando non è possibile operare una integrazione numerica.

\[ (\theta - \hat \theta) \sim Normale(\overline \theta , \space \overline \Sigma) \\ hp. \space \space \hat \theta \rightarrow \space moda \space a \space psteriori \]

Questa prima approssimazione si basa sul Teorema centrale del limite, per calcolare i due parametri della normale $\overline \theta$ e $\overline \Sigma$ basta:

$\overline \theta \rightarrow$

$\theta$

$\overline \Sigma \rightarrow$

Come secondo step per trovare $\overline \Sigma$ calcolo la derivata seconda del logaritmo della posterior e sostituisco al posto di $\theta$ la moda calcolata allo step precedente. \[ - \frac{\partial^2 \space \log ( \pi(\theta | x))}{\partial \space\theta^2} \bigg|_{\theta= \overline \theta} \]

N.B. Metodo alternativo ma molto poco utilizzato in questo ambito è calcolare $\overline \theta$ e $\overline \Sigma$ attraverso un double bootstrap.

Approssimazione di Laplace

\[hp.\space E^{\pi(.|x)}[g(\theta)]= \int_{\Theta}{g(\theta) \space \pi(\theta | x)} \space \space d\theta= \int_\Theta e^{n \space h(\theta) }\space d\theta \\ dove \space \space h(\theta)= \frac1 n \log\bigg(g(\theta) \pi(\theta | \underline x)\bigg) \]

Se faccio lo sviluppo di Taylor di $h(\theta)$ fino al secondo ordine ottengo:

\[h(\theta^*) \approx h(\theta^*) + (\theta-\theta^*)\space h^|(\theta^*)+\frac{(\theta-\theta^*)^2}{2}\space h^{||}(\theta^*) \\ dove \space \space \theta^*= \hat\theta_{SMV}= moda \space a \space posteriori \]

Quindi:

\[\int_\Theta e^{n \space h(\theta) }\space d\theta \approx \int_\Theta e^{n \space \bigg[ h(\theta^*) + (\theta-\theta^*)\space h^|(\theta^*)+\frac{(\theta-\theta^*)^2}{2}\space h^{||}(\theta^*) \bigg] }\space d\theta \approx \\ \approx e^{n \space h(\theta) } \int_\Theta e^{n \space \frac{(\theta-\theta^*)^2}{2}\space h^{||}(\theta^*) }\space d\theta \]

Tentando di riconoscere all’interno dell’integrale una Normale questa può essere riscritta come:

\[\approx e^{n \space h(\theta) } \sqrt{2 \pi} \sigma \int_\Theta \frac{1}{\sqrt{2 \pi}} \frac{1}{\sigma} e^{n \space \frac{(\theta-\theta^*)^2}{2}\space h^{||}(\theta^*) }\space d\theta \approx \\ \approx e^{n \space h(\theta) } \sqrt{2 \pi} \sigma \int_\Theta \frac{1}{\sqrt{2 \pi \sigma^2}} \space e^{-\frac12 \space \frac{(\theta-\theta^*)^2}{- \frac{1}{n\space h^{||}(\theta^*)} }\space }\space d\theta \]

Tutto cio che c’è all’interno dell’integrale equivale ad 1 in quanto è una $Normale(\theta^*, \space- \frac{1}{n\space h^{||}(\theta^*)})$, inoltre fuori dall’integrale c’è $\sigma$ che corrisponde a $\sigma=\sqrt{- \frac{1}{n\space h^{||}(\theta^*)}}$

Quindi:

\[ \approx e^{n \space h(\theta) } \sqrt{2 \pi} \sigma \approx e^{n \space h(\theta) } \sqrt{2 \pi} \sqrt{- \frac{1}{n\space h^{||}(\theta^*)}} \approx \\ \approx e^{n \space h(\theta) } \sqrt{- \frac{2 \pi}{n\space h^{||}(\theta^*)}} \]

Con questo metodo a differenza della approssimazione Normale è necessaria la conoscenza di $m(x)$ in quanto all’interno di $h$ è presente.

Metodi simulativi

I metodi di simulazione ricostruiscono le caratteristiche di una distribuzione complessa (nel nostro caso la posterior) tramite un campione di elementi indipendenti generato dalla distribuzione.

Metodi Monte Carlo

Se conosco $\pi(\theta | \underline x)$ e posso fare estazioni da essa in modo i.i.d. allora posso utilizzare le tecniche monte carlo per stimare il valore atteso della posterior. In questa parte si spiegherà il metodo utilizzato in questo ambito mentre si lascia una spiegazione dettagliata del Metodo Monte Carlo all’articolo dedicato (link articolo).

STEP 1:

$m$

$\pi(\theta | \underline x)$

\[\theta_1,...,\theta_m \sim \pi(\theta | \underline x) \\ t.c. \space \space m \gg0\]

STEP 2:

\[E^{\pi(.|x)}[g(\theta)] \cong \overline g_m = \frac{1}{m}\sum_{i=1}^m g(\theta_i) \]

N.B. Per avere una stima intervallare si può sfruttare il Teorema Centrale del Limite e la legge forte dei grandi numeri per calcolare un intervallo di confidenza:

\[ \overline g_m \space \pm \space Z_{1-\frac{\alpha}{2}} \space \space \sqrt{\hat{Var}( \overline g_m)} \\ \rightarrow \hat{Var}( \overline g_m) = \frac{\sum_{i=1}^m\bigg(g(\theta_i)- \overline g_m \bigg)}{m(m-1)}\]

Metodi Monte Carlo Importance Sampling (MCIS)

Questo secondo meccanismo di stima è una versione ponderata del metodo MC che utilizza delle funzioni ausiliarie per trovare un campione indipendente e identicamente distribuito della posterior, può essere molto utile quando non si conosce la posterior.

STEP 1:

$h(\theta)$

$\pi(\theta | \underline x)$

STEP 2:

$m$

$h(\theta)$

\[\theta'_1,...,\theta'_m \sim h(\theta) \\ t.c. \space \space m \gg0\]

$h$

STEP 3:

\[E^{\pi(.|x)}[g(\theta)] \cong \overline g_m = \frac{1}{m}\sum_{i=1}^m \frac{\pi(\theta'_i|x)}{h(\theta'_i)}g(\theta'_i) \]

Osservazione:
Se non conosco la posterior e la prior è propria allora posso intedere la posterior come $( | ) = $ da cui simulare i valori:

\[E^{\pi(.|x)}[g(\theta)]= \int_{\Theta}{g(\theta) \space \pi(\theta | x)} \space \space d\theta =\\= \int_{\Theta}{g(\theta) \space \frac{ f(\underline{x} | \theta) \pi(\theta)}{m(x)}} \space \space d\theta = \\ = \frac{ \int_\Theta g(\theta) f(\underline{x} | \theta) \pi(\theta) \space \space d\theta }{\int_\Theta f(\underline{x} | \theta) \pi(\theta) \space \space d\theta }\]

Quindi se intendo come $h(\theta)= \pi(\theta)$ (si può utilizzare anche $f(\underline x,\theta)$ se è riconoscibile), allora:

\[E^{\pi(.|x)}[g(\theta)] \cong \frac{\frac{1}{m}\sum_{i=1}^m \space \space g(\theta'_i) \space f(\underline x | \theta'_i) }{\frac{1}{m}\sum_{i=1}^m \space \space f(\underline x | \theta'_i)} = \\ = \sum_{i=1}^m\frac{ \space \space g(\theta'_i) \space f(\underline x | \theta'_i) }{ \space \space f(\underline x | \theta'_i)}\]

N.B. Anche per MCIS vale la legge forte dei grandi numeri e il teorema centrale del limite.

Monte Carlo Markov Chain (MCMC)

GIBBS SAMPLING

Impostiamo il problema classico nel caso multivariato, quindi con una serie di parametri da stimare con un approccio bayesiano:

\[x \sim f \space \space \underline \theta \space t.c. \theta_i \sim \pi_i(\theta_i) \space \space con \space \space i=1,...,k\]

Step 1: Calcolo le FULL CONDITIONAL:

$x|\theta \sim x \Rightarrow$ funzione di densità della x

$\theta_i|x \sim \pi_i(\theta_i| \underline x) \Rightarrow$

Step 2: Inizializzo un vettore in modo arbitrario:

\[(x^0, \theta_1^0,...,\theta_k^0) \Rightarrow \space in \space modo \space arbitrario\]

Step 3: Estrazione ricorsiva

Per $\theta_1$ al tempo 1

$\theta_1^1 \Rightarrow$ estraggo $\theta^1$ al tempo 1 dalla posterior $x=x^0$ e $\theta_2 = \theta_2^0$ …

Per $\theta_2$ al tempo 1

$\theta_2^1 \Rightarrow$ estraggo $\theta^2$ al tempo 1 dalla posterior $x=x^0$ e $\theta_1 = \theta_1^1, \theta_3 = \theta_3^0$ …

… così per i $\theta_i$ con $i=1,...,k$

$x^1$ estraggo $x^1$ al tempo 1 da $f(\underline x | \theta_1 = \theta_1^1, ..., \theta_k = \theta_k^1)$

Reitero il processo $m+n$ volte

Perchè per $m+n$ volte?

Perchè devo aspettare per un tempo $m$ che la catena raggiunga convergenza.

\[E^{\pi(.|x)}[g(\theta)] \cong \overline g_m = \frac{1}{n}\sum_{i=m}^{m+n} g(\theta_i) \]

METROPOLIS HASTINGS

Si tratta di un algoritmo di Accettazione e Rifiuto.

Si indica con $q(\theta,\theta')$ la proposal, ossia la legge di $\theta'$ condizionata a $\theta$

Se estraggo da $q$ un candidato $\theta'$ la sua probabilità di accettazione è:

\[\alpha(\theta, \theta')= \begin{cases} \min\bigg\{\frac{\pi(\theta'|\underline x) \space q(\theta',\theta)}{\pi(\theta|\underline x) \space q(\theta,\theta')}\space \space, 1 \bigg\} &se &\pi(\theta'|\underline x) \space q(\theta',\theta) \ne 0 \\ 1 &se &\pi(\theta'|\underline x) \space q(\theta',\theta) = 0 \end{cases}\]

Operativamente per implementare questo metodo occorre procedere seguendo 5 steps:

si fissa un $\theta^0$ in modo arbitrario (per velocizzare la converngenza si può scegliere tra media, mediano, moda …);
Si genera il candidato $\theta'$ da $q(\theta^0,\theta')$;
Si genera $u$ da una uniforme $U(0,1)$ (essa ovviamente deve essere indipendente dalle precedenti generazioni);
Se \[U \leq\alpha(\theta^0, \theta')\] Allora accetto il candidato come possibile generazione della posterior;
Reitero questo ragonamento n+m volte (per la convergenza).

Una volta ottenuto il campione casuale di generazioni della posterior opero come con i metodi precedenti, quindi:

\[E^{\pi(.|x)}[g(\theta)] \cong \overline g_m = \frac{1}{n}\sum_{i=m}^{m+n} g(\theta_i) \]

Di questo ultimo metodo ci sono due sviluppi (o modifiche) che permettono l’utilizzo di questo ache in altre situazioni:

Metropolis Random Walk

Indipendent Metropolis