Posterior synthesis

The problem we face in this article concerns Posterior synthesis. Let us state the problem in the most general way possible:

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

The main methods that can be used to solve the problem just stated can be divided into three categories:

Numerical integration methods

Analytical methods

Simulation methods

Numerical Integration Methods

Numerical integration consists of a series of methods that estimate the value of a definite integral without having to calculate the antiderivative of the integrand function.

The need for numerical integration arises because not all functions admit an antiderivative in explicit form (for example the Gaussian curve), or the antiderivative can be very complicated to evaluate, or we are dealing with a function available only at certain points.

Numerical integration methods can be distinguished into two macro-categories:

Newton-Cotes formulae (e.g. trapezoidal method);
Gauss formulae (e.g. Monte Carlo method).

These methods for posterior approximation work well only when the parameter space is very small; in this context they are not used much, although for conjugate posteriors Monte Carlo methods are an excellent choice.

Furthermore, often we do not know $m(x)= \int_{\Theta}{f(\underline{x} | \theta) \pi(\theta)} \space \space d\theta$ and this forces us to use approximations.

Analytical methods

Normal Approximation

The Normal approximation is a method that uses the asymptotic properties of the Normal distribution to solve the Posterior synthesis problem when numerical integration is not feasible.

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

This first approximation is based on the Central Limit Theorem. To calculate the two parameters of the normal $\overline \theta$ and $\overline \Sigma$ it suffices to:

$\overline \theta \rightarrow$

$\theta$

$\overline \Sigma \rightarrow$

As the second step, to find $\overline \Sigma$ calculate the second derivative of the logarithm of the posterior and substitute in place of $\theta$ the mode calculated in the previous step. \[ - \frac{\partial^2 \space \log ( \pi(\theta | x))}{\partial \space\theta^2} \bigg|_{\theta= \overline \theta} \]

N.B. An alternative method, though rarely used in this context, is to calculate $\overline \theta$ and $\overline \Sigma$ through a double bootstrap.

Laplace Approximation

\[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 \\ where \space \space h(\theta)= \frac1 n \log\bigg(g(\theta) \pi(\theta | \underline x)\bigg) \]

If we perform the Taylor expansion of $h(\theta)$ up to the second order we get:

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

Therefore:

\[\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 \]

Trying to recognise a Normal distribution within the integral, this can be rewritten as:

\[\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 \]

Everything inside the integral equals 1 since it is a $Normale(\theta^*, \space- \frac{1}{n\space h^{||}(\theta^*)})$, furthermore fuori dall’integrale c’è $\sigma$ which corresponds to $\sigma=\sqrt{- \frac{1}{n\space h^{||}(\theta^*)}}$

Therefore:

\[ \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^*)}} \]

With this method, unlike the Normal approximation, knowledge of $m(x)$ in quanto all’interno di $h$ it is present.

Simulation methods

Simulation methods reconstruct the characteristics of a complex distribution (in our case the posterior) through a sample of independent elements generated from the distribution.

Monte Carlo Methods

If we know $\pi(\theta | \underline x)$ and posso fare estazioni da essa in modo i.i.d. allora posso utilizzare the tecniche monte carlo per estimatere the value atteso of the posterior. In this parte it spiegherà the method utilizzato in this ambito mentre it lascia a spiegazione dettagliata of the 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. To obtain an interval estimate, we can exploit the Central Limit Theorem and the Strong Law of Large Numbers to calculate a confidence interval:

\[ \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)}\]

Monte Carlo Methods Importance Sampling (MCIS)

This second estimation mechanism is a weighted version of the MC method that uses auxiliary functions to find an independent and identically distributed sample of the posterior; it can be very useful when the posterior is unknown.

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) \]

Observation:
If we do not know the posterior and the prior is proper, then we can express the posterior as $( | ) = $ from which to simulate the values:

\[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 }\]

Therefore if we take as $h(\theta)= \pi(\theta)$ (one can also use $f(\underline x,\theta)$ if it is recognisable), then:

\[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. For MCIS too, the Strong Law of Large Numbers and the Central Limit Theorem apply.

Monte Carlo Markov Chain (MCMC)

GIBBS SAMPLING

Let us set up the classical problem in the multivariate case, with a series of parameters to estimate using a Bayesian approach:

\[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: Compute the FULL CONDITIONAL:

$x|\theta \sim x \Rightarrow$ density function of x

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

Step 2: Initialise a vector arbitrarily:

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

Step 3: Recursive extraction

For $\theta_1$ at time 1

$\theta_1^1 \Rightarrow$ I draw $\theta^1$ at time 1 from the posterior $x=x^0$ and $\theta_2 = \theta_2^0$ …

For $\theta_2$ at time 1

$\theta_2^1 \Rightarrow$ I draw $\theta^2$ at time 1 from the posterior $x=x^0$ and $\theta_1 = \theta_1^1, \theta_3 = \theta_3^0$ …

… and so on for the $\theta_i$ con $i=1,...,k$

$x^1$ I draw $x^1$ at time 1 da $f(\underline x | \theta_1 = \theta_1^1, ..., \theta_k = \theta_k^1)$

I reiterate the process $m+n$ times

Why for $m+n$ times?

Because we must wait for a time $m$ for the chain to reach convergence.

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

METROPOLIS-HASTINGS

This is an Accept-Reject algorithm.

We denote by $q(\theta,\theta')$ the proposal, i.e. the law of $\theta'$ conditioned on $\theta$

Se I draw da $q$ a candidate $\theta'$ its acceptance probability is:

\[\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}\]

Operationally, to implement this method we proceed following 5 steps:

fix a $\theta^0$ arbitrarily (to speed up convergence one can choose among mean, median, mode...);
Generate the candidate $\theta'$ da $q(\theta^0,\theta')$;
Generate $u$ from a uniform $U(0,1)$ (it must obviously be independent from previous generations);
Se \[U \leq\alpha(\theta^0, \theta')\] Then accept the candidate as a possible generation from the posterior;
Reitero this ragonamento n+m times (per the convergenza).

Once the random sample of posterior generations is obtained, we proceed as with the previous methods:

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

Di this ultimo method ci are due sviluppi (o modifiche) that permettono l’utilizzo di this ache in altre situazioni:

Metropolis Random Walk

Independent Metropolis