title: "Approximations"

author: "Adrien Osakwe"

format: html

editor: visual

Previously we have looked at Bayesian models with closed form solutions for the marginal likelihood and posterior. However, in cases where more complex models are required, the marginal likelihood is often not closed form making it impossible to solve the posterior analytically.

In this cases, we are required to **approximate** the posterior from which we can then compute the statistics of interest which we saw in the previous section.

There are two key ways to approximate the posterior:

1. Sampling/Simulation: we can generate random samples from the posterior distribution to then use the empirical results as an estimate of the posterior.

2. Variational Inference: we can approximate the posterior by using a closed-form surrogate distribution.

A simple form of variational inference is normal approximation, where the central limit theorem is used to justify the use of the normal distribution as a surrogate posterior. In this chapter we will focus on sampling methods for posterior approximation.

To approximate the posterior by simulation, we need to generate random samples from the posterior. This is trivial if the posterior is a known distribution with a derived simulation method. However, this becomes interesting when the posterior **isn't** closed form. We can achieve this by making use of the unnormalized posterior density or any function $q(\theta;y)$ that is proportional to it.:

and generating random samples from this to approximate the posterior.

To generate the random samples, we can use **rejection sampling** or **importance sampling** for simple models. There are also many probabilistic programming tools available to automatically generate simulations.

This method, also know as **direct discrete approximation**, works as follows:

1. Create an even-spaced grid $g_1 = a + i/2,…,g_m = b - i/2$ where $a$ and $b$ are the lower and upper bounds of the parameter space of the posterior and $i$ and $m$ are the grid increments and number of grid points respectively.

2. Evaluate the values of the unnormalized posterior density across the grid points $q(g_j;y)$ and normalize them by their sum to estimate the posterior values:

$$

3. For every iteration $s = 1,…,S:$

1. Generate a sample for $\theta_s$ from a categorical distribution with $m$ outcomes (\$g_1,...,g_m\$) with respective probabilities $\hat{p}_1,...,\hat{p}_m$

2. Add zero-centered uniformally distributed noise $\epsilon$ that has an interval length equal to the grid interval spacing such that:

$$\hat{\theta}_s = \theta_s + \epsilon$$

where $\epsilon \sim Uniform(-i/2,i/2)$

This generates a process similar to numerical integration, meaning that this approximation approach is limited to a finite interval.

Note: I have been strictly using $\theta$ as the symbol for parameters across the different models used to make it clear what symbol represents the distribution parameter. In practice, different distributions make use of different symbols for their parameter. Examples are $\alpha, \beta, \lambda$ some of which were used in previous sections when defining the prior distributions.

We can demonstrate the usage of grid approximation with the Poisson-gamma model we explored previously. Simulation isn't necessary here since the posterior is closed-form, but this will at least demonstrate how grid approximation is able to approximate it.

Recall that the posterior of the Poisson-Gamma model had the following form:

Note that we needed to explore $\theta$ on a fixed interval (0,25) to be able to approximate it. We can determine the plausible range of values to explore using the prior (if it is informative) or the likelihood. The key requirement is that we make sure the interval covers nearly all of the pmf for the distribution.

We can clearly see that the grid approximation has allowed us to correctly approximate the posterior without having to evaluate it analytically.

Let's look at an example of a non-conjugate distribution, where the normalizing constant is such that the posterior is not a known distribution and is therefore **not closed-form (intractable)**. A common example is using a log-normal prior for a Poisson likelihood. In essence, if our random variable is $X$ and is normally distributed, we can use a reparameterization of the r.v, $Y = e^X$ to get a log-normal distribution with the **mean and variance being equal to the log of the mean and variance for** $X$ .

Our model therefore uses the following distributions:

The pdf for the prior is as follows:

which gives the following un normalized posterior pdf

And a marginalized joint pdf (numerator) we cannot integrate over (intractable). We can now apply grid approximation! We will be reusing the distribution we create in the last example.

We can see that our approximation is close to the ground truth, but is slightly off. This could be seen as an indication that the selected prior may not be the best option (see result with gamma prior above)!

We have now seen that a sufficient number of simulations can approximate the ground truth posterior. These simulations can also be used to calculate summary statistics such as the posterior mean, variance and credible intervals!

In general, computing integrals via simulations is referred to as **Monte Carlo integration** or the **Monte Carlo method**. These types of methods depend on the **strong law of large numbers (SLL)**.

Where we have a sequence of random variables $Y_i,…,Y_n$ that are i.i.d. with a finite expected value $\mu$ , we have that almost surely (a.s.):

We are able to say almost surely as the above statement is said to converges to the result with a probability of one.

SLL implies that the sample mean of an i.i.d sequence of random variables will always converge to the expected value of the underlying distribution with sufficient sample size. In the case of coin tosses, we would therefore expect that after enough trials, by SLL, that the expected value is 0.5 (assuming it is an unbiased coin!).

We can revisit our example from the beginning of the chapter where we estimate $\theta$. By SLL, we expect that that average of the simulated values will converge to the true posterior mean of $\theta$ provided we simulate enough values. We can therefore approximate the posterior expectation using this mean. In this case, we know the posterior expectation of the Poisson-gamma model

so we can validate our approximation.

The same can be done to approximate the true posterior variance.

We can also approximate posterior probabilities within an interval defined by an indicator $I_{(a,b)}(x)$ where $I_{(a,b)}(x) = 1$ if $x \in (a.b)$ and $I_{(a,b)}(x) = 0$ otherwise:

Therefore, we are also able to estimate quantiles:

In cases where we can't generate samples from the target distribution, Monte Carlo methods allow us to generate viable samples by means of a surrogate distribution. In essence, we aim to identify a pdf $f_0$ for which sampling is possible to then estimate the features of interest of our target distribution's pdf $f$. For this to work, it is imperative that **the** **selected pdf has the support for the same domain as the target pdf**. In the case of determining an expected value we get:

Which gives us

Using the SLL as discussed previously, we can then use the **estimator** $\hat{I}^{(f_0)}_N(g)$ to get the expected value over the target pdf:

In practice, $\hat{I}^{(f_0)}_N(g)$ is known as the **importance sampling** estimator. To emphasize the notion of importance, we can rearrange the equation as follows:

Where $w_0(X_i)$ represents the **importance sampling weight**. We would choose to use the importance sampler estimator over the standard estimator in cases where its variance is less than that of the standard estimator.

Let's look at an example. Let us say that we want to get the following expectation over $f(x)$

Where we know that $f(x)$ is the pdf of a Gamma distribution $X \sim Gamma(2,3)$. We aim to compute $\mathbb{E}[X^3]$. As we know the ground truth, we know that the expected value for the random variable of a gamma function is:

We can now evaluate importance sampling estimators that use a series of $f_0$ functions:

- $Gamma(2,3)$ - Standard MC Estimator

- $Gamma(6,3)$

- $Student(5)$ centered around $x = 2$

- $Normal(3,2)$

- $Gamma(2,6)$

Doing this for multiple runs will help us determine the variance of each estimator.

As we can see, the standard MC estimator actually has a lot more uncertainty than the IS approaches. We notice that the estimators using a gamma pdf work quite well.

When choosing $f_0(x)$, we want to ensure that the variance is finite so that our confidence on the estimate is quantifiable. The estimator's variance is finite if and only if

has finite variance. In other words, it is finite if

is finite. We therefore conclude that an **optimal** choice for $f_0$ occurs when$f_0(x) \propto |g(x)|f(x)$

When the ratio $\frac{f(x)}{f_0(x)}$ has unbounded support, the variance of our estimator will not be finite. We therefore need a way to ensure that the ratio remains bounded at the tails for $f$ with unbounded support. In practice, we find that

allow us to derive the following estimator

Which, if it exists, may have a smaller variance than the variance we saw earlier.

-Note: consider adding something about harmonic mean derivation

Despite being similar, there is a key difference between importance and rejection sampling:

1. Importance sampling focuses on approximating numerical integration with respect to $f(x)$

2. Rejection sampling focuses on generating i.i.d samples from $f(x)$

Based on these differences, we can consider **Sampling Importance Resampling** (SIR) which can sample from density $f(x)$ by *re-weighting* and then *resampling* samples from $f_0(x)$:

1. Generate samples $x_1,…,x_n$ from $f_0$

2. Calculate re-normalized weights $w_1,…,w_n$ $$w_i = \frac{\frac{f(X_i)}{f_0(X_i)}}{\sum_{j = 1}^N\frac{f(X_j)}{f_0(X_j)}} $$

3. Generate new samples (resample) from the discrete distribution of the samples $x$ using the weights $w$ as the probability masses.

Note: make an example

Grid approximation worked well in our 1-dimensional example. However, this is not always the case when we begin to work with high-dimensional data. To cover the parameter space we defined in our example, we generated a grid with 2000 points. However, if we had additional dimensions each with 2000 points, we would find ourselves working with $2000^d$ grid points. This becomes quite impractical in what has been a very simple example. Clearly, grid approximation will not be viable for higher dimensions.

As a result, most sampling approaches for more complex models usually make use of **Markov Chain Monte Carlo (MCMC)** methods. This family of methods focuses on iterative sampling from a markov chain for which the distribution at convergence matches the target distribution. In our case, this is the posterior distribution.

A Markov chain is a sequence of random variables $X_1,X_2,…$ that have a Markov property. That is, for any discrete time point $i$ :

That is to say, for any discrete time point \$i\$, the state of the random variable at the next time point, $X_{i+1}$ depends **only** on the current state $X_i$ . The set of all possible values or states that the random variables can take is known as the **state space** $S$.

Simplistic sampling methods such as the ones we saw previously generate i.i.d. samples from the target distribution. However, the values we sample from the state space S ($\theta_1,…,\theta_S$) by an MCMC method experience high auto-correlation. That is, the value sampled at the next time point will be similar to the current value. Although this would be an issue for a small sample size, MCMC overcomes this by generating a large number of samples and using the full sample to approximate the posterior distribution. MCMC algorithms converge when they reach the stationary distribution. That is, the stationary distribution is a distribution $\pi(x)$ achieved at a state where for all $i \in \mathbb{Z}^+$ :

which implies that the next state of the random variable is the same as the current state. Once it has reached convergence, sampling from the stationary distribution will eventually lead to an approximation of the true posterior. Determining when convergence has been reached is difficult to determine. There are different diagnostics available to explore this, one of which is to running different MCMC instances with different initializations and see if they reach a similar distribution.

As the first iterations of MCMC sampling do not represent the posterior yet, they are usually discarded. The number of iterations that are discarded is up to the user and is referred to as the **burn-in** period. Markov chains that were designed to estimate target posterior distributions are known as **MCMC samplers**. The two most popular examples are the **Gibbs sampler** and the **Metropolis-Hastings sampler**. The Gibbs sampler is actually a special case of the MH sampler.

The Gibbs sampler works by incrementally updating each of the components of the parameter vector $\theta$. We can assume that the parameter vector is multi-dimensional where for each component $\theta_j$ , the sampler generates a value using the conditional posterior distribution of the component given all the other components.

Notice that the current component **is omitted** from the parameter vector being used as prior information!

We can explore the concept with a two-dimensional example. We assume that we have one observation $(y_1,y_2) = (0,0)$ which comes from a two-dimensional normal distribution $N(\mu,\Sigma_0)$ , where the parameter we are interested in is the mean vector $\mu = (\mu_1,\mu_2)$ and the covariance matrix

$$\Sigma_0 =

\mu_1|\mu_2,Y \sim N(y_1 + \rho(\mu_2-y_2),1 - \rho^2)

\mu_2|\mu_1,Y \sim N(y_2 + \rho(\mu_1-y_1),1 - \rho^2)

A = min(1,\frac{\pi(z)Q(x|z)}{\pi(x)Q(z|x)}

p(y'|y) = \int p(y'|\theta)p(\theta|y)d\theta