# Posterior Approximation

## The Law of Large Numbers and Monte Carlo Estimation

What if there's no well-defined PDF for our product of likehood and prior(s)? It means we cannot use `scipy.stats.norm`, and its convenient `mean()`, `pdf()`, `cdf()`, or `ppf()` functions to establish probabilistic answers to questions about the posterior. Is it possible to answer such questions in a different way? Assuming we have a sample $y_1,\ldots,y_n$ from an _unknown_ distribution:


```python
import numpy as np

y = np.random.normal(loc=130, scale=10, size=10)  # although using samples from a normal distribution, pretend to not know the true distribution
```

The distribution is not defined by a PDF, but by a **sample**. According to the [Central Limit Theorem](https://en.wikipedia.org/wiki/Central_limit_theorem), it is possible to approximate the distribution mean with the sample mean, with large enough sample. Is $n=10$ good enough?


```python
np.mean(y)
```


```python
np.std(y)
```

As [mentioned before](./notebooks/01_bayes_rule_intro.ipynb#Continuous:), the probability of $Y$ being between e.g. 110 and 130 is

$$Pr(110 \le Y \le 130) = \int_{110}^{130}p(y)$$

However, $p(y)$ is undefined and the only thing we have is a sample. But the sample can be used to approximate the integral:


```python
def sample_cdf(x, sample):
    return (sample <= x).sum() / sample.size  # "integrating" by counting
```


```python
sample_cdf(130, y) - sample_cdf(110, y)
```

This can be compared with the probability obtained from the true normal distribution:


```python
from scipy.stats import norm

unknown_norm = norm(loc=130, scale=10)
unknown_norm.cdf(130) - unknown_norm.cdf(110)
```

Samples also allow us to compute quantiles:


```python
def sample_ppf(p, sample):
    p_index = int(np.round((sample.size - 1) * p))
    return np.sort(sample)[p_index]
```

So the range of $y$ for which the probability is 80% that it contains $y$'s true mean is:


```python
f'[{sample_ppf(0.1, y)} - {sample_ppf(0.9, y)}]'
```

where the interval from the true normal distribution would be:


```python
f'[{unknown_norm.ppf(0.1)} - {unknown_norm.ppf(0.9)}]'
```

Using a larger sample leads to more accurate estimations:


```python
y_large = np.random.normal(loc=130, scale=10, size=10000)

{
    'mean (130)': np.mean(y_large),
    'std (10)': np.std(y_large),
    'p_110_130 (0.4772)': sample_cdf(130, y_large) - sample_cdf(110, y_large),
    '80% ([117.2 - 142.8])': f'[{sample_ppf(0.1, y_large)} - {sample_ppf(0.9, y_large)}]'
}


```

### Chaining Simulations for Hierarchical Models

This process of computing means, variances, probabilities, percentiles (or any other quantities of interest) from simulated samples is called **Monte Carlo Estimation**. When using more complex models, simulations such as above can also be **chained**. For example, consider a model similar to the one in the [previous note](./02_generative_models.ipynb#Example,-observing-a-single-value-with-a-simple-model):

<div style="font-size: 2em">
$$
\begin{align}
\mu &\sim\color{red}{\mathcal{t}(m_0, s_0, \nu)}\,\mathrm{(prior)}\\
y &\sim\color{blue}{\mathcal{N}(\mu, \sigma^2_0)}\,\mathrm{(likelihood)}\\
m_0 &=130\\
s_0 &=10\\
\nu &=(n-1)\,\text{degrees of freedom for a sample size of}\,n\\
\sigma^2_0 &=100
\end{align}
$$
</div>

The joint density $p(y, \mu) = \color{blue}{p(y\mid\mu)}\color{red}{p(\mu)}$ cannot be expressed as a known family (e.g. Normal) of probability densities, but can be estimated by simulation:


```python
from scipy.stats import t

# basing the prior on a sample of n=10, simulating m=1000 samples
mu = t.rvs(df=9, loc=130, scale=10, size=1000)

# an array of n location parameters can be plugged in to get n new samples
y_chained = norm.rvs(loc=mu, scale=10)

{
    'mean (130)': np.mean(y_chained),
    'std (10)': np.std(y_chained),
    'p_110_130 (0.4772)': sample_cdf(130, y_chained) - sample_cdf(110, y_chained),
    '80% ([117.2 - 142.8])': f'[{sample_ppf(0.1, y_chained)} - {sample_ppf(0.9, y_chained)}]'
}
```

These are samples from the **prior predictive** distribution, i.e. given the likelihood and all the priors in a model, what values of $y$ can we expect (before having observed any data)?

- Can we use the above approach to sample from the _posterior_ distribution $p(\mu\mid y)$? Why (not)?

## Markov Chains

A [**Markov Chain**](https://en.wikipedia.org/wiki/Markov_chain) is a sequence of random variables $X_1,\ldots,X_n$, with $1,\ldots,n$ representing points in time, where the probability (density) of $X_{t+1}$ only depends on the value $X_t$, i.e.

$$
p(X_{t+1}\mid X_t,X_{t-1},\ldots,X_2,X_1) = p(X_{t+1}\mid X_t)
$$

A random walk with $X_0\sim\mathcal{N}(0, 1)$ and $X_{t+1}\mid X_t\sim\mathcal{N}(X_t, 1)$ is an example of a Markov Chain:


```python
def normal_random_walk(sigma, mu_init=0):
    mu = mu_init
    
    while True:
        yield mu
        mu = np.random.normal(loc=mu, scale=sigma)
```


```python
import plotly.graph_objs as go

random_walk_iter = normal_random_walk(1)
n=1000

random_walk_trace = go.FigureWidget(
    data=[go.Scatter(x=list(range(n)), y=[next(random_walk_iter) for _ in range(n)], showlegend=False)]
)

random_walk_trace
```

What happens with multiple random walks?


```python
import pandas as pd

trace_df = pd.DataFrame()

for i in range(10):
    random_walk_iter = normal_random_walk(1)  # re-initialize generator for every new random walk
    trace = [next(random_walk_iter) for _ in range(n)]
    trace_df[f'trace_{i}'] = trace
    random_walk_trace.add_scatter(x=list(range(n)), y=trace, showlegend=False)
```

They look ... quite random:


```python
trace_df.describe()
```

However, with a little tweak of the random walk...


```python
def tweaked_random_walk(sigma, mu_init=0, phi=0.9):
    mu = mu_init
    
    while True:
        yield mu
        mu = np.random.normal(loc=phi*mu, scale=sigma)
```


```python
n=1000
sigma=20

random_walk_iter = tweaked_random_walk(1, mu_init=np.random.normal(loc=0, scale=sigma))


random_walk_trace = go.FigureWidget(
    data=[go.Scatter(x=list(range(n)), y=[next(random_walk_iter) for _ in range(n)], showlegend=False)]
)

random_walk_trace
```


```python
trace_df = pd.DataFrame()

for i in range(10):
    random_walk_iter = tweaked_random_walk(1, mu_init=np.random.normal(loc=0, scale=sigma))  # re-initialize generator for every new random walk
    trace = [next(random_walk_iter) for _ in range(n)]
    trace_df[f'trace_{i}'] = trace
    random_walk_trace.add_scatter(x=list(range(n)), y=trace, showlegend=False)
    
trace_df.describe()
```

... the traces converge to a **stationary distribution** (with mean 0). When experimenting with different values for $\phi$ in the `tweaked_random_walk()` generator, it turns out that the random walks converge in the case of $-1\lt\phi\lt1$. In that case, the generator `tweaked_random_walk()` can be considered a **sampler** for $\mathcal{N}(0, \frac{1}{1-\phi^2})$.

**Note:** Markov Chains are not limited to continuous variables. It is possible to create discrete Markov Chains by modeling probabilistic state transitions between discrete variables.

Would it be possible to construct such a sampler for arbitrary complex distributions? Such as a posterior $p(\theta\mid y)$?

## Getting Picky - The Metropolis Sampler

Assume a model that estimates the month-over-month _change_ in mean exercise heart rate $\mu$ for a population. Such a model can be defined as:

<div style="font-size: 2em">
$$
\begin{align}
y &\sim\color{blue}{\mathcal{N}(\mu, 1)}\,\mathrm{(likelihood)}\\
\mu &\sim\color{red}{\mathcal{t}(0, 1, 1)}\,\mathrm{(prior)}\\
\end{align}
$$
</div>

with the prior for $\mu$ coming from a standard $\mathcal{t}$ distribution, i.e. centered around 0 with scale equal to 1 standard deviation. This allows to write the posterior as:

<div style="font-size: 2em">
$$
\begin{align}
p(\mu\mid y_1,\ldots,y_n) &\propto \color{blue}{p(y_1,\ldots,y_n\mid\mu)}\color{red}{p(\mu)}\\
&= \color{blue}{\prod_{i=1}^n \frac{1}{\sqrt{2\pi}}exp\left\{-\frac{1}{2}{(\color{black}{y_i} - \mu)}^2\right\}}\color{red}{\frac{1}{\pi(1+\mu^2)}}\\
&\propto\ldots\\
&\propto \frac{exp\left[n(\bar{y}\mu - \frac{\mu^2}{2})\right]}{1 + \mu^2}\,= g(\mu)
\end{align}
$$
</div>

where $n$ and $\bar{y}$ are known and fixed. To summarize, $g(\mu)$ is:

- a PDF,
- proportional to $p(\mu\mid y)$

The [Metropolis algorithm](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm) allows to sample from the posterior $p(\mu\mid y)$ by creating a Markov Chain with a stationary distribution that approximates the posterior. It does so by iteratively accepting or rejecting samples drawn from another distribution that is easier to sample. 

Remember this?


```python
norm.pdf(2) / norm.pdf(1)
```

<div style="font-size: 2em">
$$
\begin{align}
p(\mu\mid y) &\propto g(\mu) \Leftrightarrow\\
\frac{p(\mu_1\mid y)}{p(\mu_2\mid y)} &= \frac{g(\mu_1)}{g(\mu_2)}\\
\text{or}\\
\frac{p(\mu_{t+1}\mid y)}{p(\mu_t\mid y)} &= \frac{g(\mu_{t+1})}{g(\mu_t)}\\
\end{align}
$$
</div>

Since we can evaluate $g(\mu)$, the ratio $\frac{g(\mu_1)}{g(\mu_2)}$ tells us how much more/less likely it is to see $\mu_1$ compared to $\mu_2$ in the **target distribution** $p(\mu\mid y)$! In the Metropolis sampler, this ratio is used to accept or reject new samples $\mu_{t+1}$ based on the previous sample $\mu_t$. This approach can be integrated into the random-walk sampler:


```python
def log_g(y_bar, n, mu):
    """
    Returns the joint density of y and mu assuming a normal likelihood and t prior,
    at log scale (to prevent under/overflow issues)
    """
    return n * (y_bar * mu - mu**2 / 2) - np.log(1 + mu**2)
```


```python
def metropolis_sampler(y, mu_init=0.0, sigma=1):
    samples = 0  # numer of samples drawn so far
    accept = 0  # number of accepted samples so far
    mu = mu_init  # current state of sampler
    y_bar = np.mean(y)
    n = y.size
    
    while True:
        yield (mu, accept/samples if samples > 0 else 0)
        
        candidate_mu = np.random.normal(loc=mu, scale=sigma)
        acceptance_ratio = np.exp(log_g(y_bar, n, candidate_mu) - log_g(y_bar, n, mu))
        
        if np.random.uniform() < acceptance_ratio:
            # accept the new candidate
            mu = candidate_mu
            accept += 1
        
        samples += 1
    
```

Now, assume we observe the following month-over-month changes of mean exercise heart rate for 10 athletes:


```python
y = np.array([3, 0, -2, 5, 1, 0, 8, -4, -1, 3])
np.mean(y)
```


```python
sampler = metropolis_sampler(y, sigma=5)
samples = 1000

large_step_trace = [next(sampler) for _ in range(samples)]
np.mean([sample[0] for sample in large_step_trace])
```


```python
metropolis_trace = go.FigureWidget(
    data=[go.Scatter(x=list(range(samples)), y=[sample[0] for sample in large_step_trace], name=f'large steps, acceptance rate = {large_step_trace[-1][1]:.3f}')]
)

metropolis_trace
```


```python
sampler = metropolis_sampler(y, sigma=0.05)
small_step_trace = [next(sampler) for _ in range(samples)]
metropolis_trace.add_scatter(x=list(range(samples)), y=[sample[0] for sample in small_step_trace], name=f'small steps, acceptance rate = {small_step_trace[-1][1]:.3f}');
np.mean([sample[0] for sample in small_step_trace])
```


```python
sampler = metropolis_sampler(y, sigma=1.2)
good_step_trace = [next(sampler) for _ in range(samples)]
metropolis_trace.add_scatter(x=list(range(samples)), y=[sample[0] for sample in good_step_trace], name=f'good step size, acceptance rate = {good_step_trace[-1][1]:.3f}');
np.mean([sample[0] for sample in good_step_trace])
```

- Why is the posterior mean for $\mu$ smaller than the sample mean $\bar{y}$?

Finally, let's take a look at the posterior distribution for $\mu$:


```python
go.FigureWidget(
    data=[go.Histogram(x=[sample[0] for sample in good_step_trace])]
)
```

3d0214c04b87fe5f35c63cb55cb6a0601e42490d
ipynb
Jupyter Notebook
03_sampling_markov_chains.ipynb
jsamoocha/bayes_notes
cb1ba474e687e1d34ca2e9806ce1492c4b15b34f
MIT
2020-02-12T02:29:59.000Z
2020-02-12T02:29:59.000Z
03_sampling_markov_chains.ipynb
jsamoocha/bayes_notes
cb1ba474e687e1d34ca2e9806ce1492c4b15b34f
MIT
2020-03-24T16:26:40.000Z
2021-06-01T22:58:32.000Z
03_sampling_markov_chains.ipynb
jsamoocha/bayes_notes
cb1ba474e687e1d34ca2e9806ce1492c4b15b34f
MIT
Qwen/Qwen-72B
1. YES
2. YES