$\newcommand{\pr}{\textrm{Pr}}$
$\newcommand{\l}{\left}$
$\newcommand{\r}{\right}$
$\newcommand\given[1][]{\:#1\vert\:}$
$\newcommand{\var}{\textrm{Var}}$
$\newcommand{\mc}{\mathcal}$
$\newcommand{\lp}{\left(}$
$\newcommand{\rp}{\right)}$
$\newcommand{\lb}{\left\{}$
$\newcommand{\rb}{\right\}}$
$\newcommand{\iid}{\textrm{i.i.d. }}$
$\newcommand{\ev}{\textrm{E}}$
$\newcommand{\odds}{\textrm{odds}}$

So far, we've looked at simple posterior distributions that we can evaluate
exactly. However, sometimes we may be interested in more complex posteriors
that we cannot evaluate numerically. We may also be interested in other quantities
like the $\pr\lp \theta \in A \given y_1, \cdots, y_n\rp$ for arbitrary sets $A$.
We may also want to calculate the posterior of some function of $\theta$ 
like $\l|\theta_1 - \theta_2\r|$. We can use Monte Carlo sampling to do this.

# 4.1 The Monte Carlo method

**Monte Carlo approximation** is based on random sampling. Let $\theta$ be a parameter
of interest and let $y_1, \cdots, y_n$ be the numerical values of a sample from a distribution
$p\lp y_1, \cdots, y_n \given \theta \rp$. Suppose we take $S$ independent, random samples
from the posterior distribution $p\lp \theta \given y_1, \cdots, y_n \rp$:
$$\theta^{\lp 1\rp}, \theta^{\lp 2\rp}, \cdots, \theta^{\lp S\rp} \sim \iid p\lp \theta \given y_1, \cdots, y_n \rp$$
Since we are sampling from the posterior, these $S$ samples will approximate the posterior,
and the accuracy will increase as $S$ increases. The empirical distribution 
$\lb\theta^{\lp 1\rp}, \theta^{\lp 1\rp}, \cdots, \theta^{\lp S\rp}\rb$ is known as the 
**Monte Carlo approximation** to $p\lp \theta \given y_1, \cdots, y_n \rp$. As $S\rightarrow\infty$,

* $\bar{\theta} = \sum_{s=1}^S \theta^{\lp s\rp} / S \rightarrow \ev\l[\theta \given y_1, \cdots, y_n\r]$ (the empirical mean of the samples goes to the expected value of the posterior)
* $\sum_{s=1}^S \lp\theta^{\lp s\rp} - \bar{\theta}\rp^2 / \lp S-1\rp \rightarrow \var\l[\theta \given y_1, \cdots, y_n\r]$ (the empirical variance of the sample goes to the variance of the posterior)
* #$\lp \theta^{\lp S\rp} \leq c\rp / S \rightarrow \pr\lp\theta \leq c \given y_1, \cdots, y_n\rp$ (the empirical
probability of observing a sample less than $c$ approximates the cdf of the posterior)
* the empirical distribution of $\lb \theta^{\lp 1\rp}, \cdots, \theta^{\lp S\rp} \rb
\rightarrow p\lp \theta \given y_1, \cdots, y_n \rp$
* the median of $\lb \theta^{\lp 1\rp}, \cdots, \theta^{\lp S\rp} \rb \rightarrow \theta_{1/2}$
* the $\alpha$-percentile of $\lb \theta^{\lp 1\rp}, \cdots, \theta^{\lp S\rp} \rb \rightarrow \theta_\alpha$

Almost any aspect of a posterior distribution can be approximated arbitrarily exactly using Monte Carlo,
and this is true for almost all functions $g\lp\theta\rp$ as well.

You can look at plots of estimated values to see when you have performed enough Monte Carlo 
samples. For instance, you could plot the empirical mean of the samples (which would be the
expected value of the posterior distribution). With enough samples, this estimate should
converge to some value.

If you are using the mean of the Monte Carlo samples $\bar{\theta}$ to estimate the mean of the posterior,
the Central Limit Theorem says that $\hat{\theta}$ is approximately normally distributed with mean
$\ev\l[\theta\given y_1,\cdots,y_n\r]$ and standard deviation $\sqrt{\var\l[\theta\given y_1,\cdots,y_n\r] / S}$.
The Monte Carlo standard error is the approximation to this standard deviation $\sqrt{\hat{\sigma}^2 / S}$ where
$\hat{\sigma}$ is the empirical standard deviation. An approximate 95% confidence interval for the posterior
mean of $\theta$ is $\hat{\theta} \pm 2\sqrt{\hat{\sigma}^2/S}$. You can choose $S$ large enough to get the
desired precision for $\theta$.

# 4.2 Posterior inference for arbitrary functions

Suppose we are interested in the posterior distribution of some computable function $g\lp\theta\rp$ 
of $\theta$. For the binomial model, we may be interested in the log odds
$$\log \odds\lp\theta\rp = \log\frac{\theta}{1-\theta} = \gamma$$
The law of large numbers says that if we sample from the posterior distribution of $\theta$,
then the average value of $\log\frac{\theta^{\lp S\rp}}{1-\theta^{\lp S\rp}}$ converges to
$\ev\l[\log\frac{\theta}{1-\theta}\given y_1, \cdots, y_n\r]$. You can also calculate
$g\lp\theta^{\lp 1\rp}\rp, \cdots, g\lp\theta^{\lp S\rp}\rp$ independently as you sample.
From these samples, as $S\rightarrow\infty$, you get all of the information from the bullet
points above (mean, variance, empirical distribution, median, $\alpha$-percentile).

# 4.3 Sampling from predictive distributions

The predictive distribution of a random variable $\tilde{Y}$ is a probability distribution for
$\tilde{Y}$ such that 

* known quantities have been conditioned out
* unknown quantities have been integrated out

A predictive distribution that integrates over unknown parameters but is not conditional
on observed data is called a **prior predictive distribution**. Such a distribution
can be useful for evaluating whether a prior distribution for $\theta$ actually
translates into reasonable prior beliefs for observable data $\tilde{Y}$.

After we observe data, the **posterior predictive distribution** is
$$\pr\lp\tilde{Y} = \tilde{y} \given Y_1=y_1, \cdots, Y_n=y_n\rp = 
\int p\lp\tilde{y}\given\theta\rp p\lp\theta\given y_1, \cdots, y_n\rp d\theta$$
In the case of a Poisson model with a gamma prior, the posterior predictive distribution
was a negative binomial.

Often, we can sample from $p\lp\theta\given y_1, \cdots, y_n\rp$ and $p\lp y\given\theta\rp$
but it is too complicated to sample from the posterior predictive distribution. We can use
a Monte Carlo procedure to sample indirectly from the posterior predictive distribution instead.
So we want a set of samples from $\tilde{Y}$ from its posterior  predictive distribution. We can
do this by sampling 
\begin{align}
&\textrm{sample } \theta^\lp 1\rp \sim p\lp\theta\given y_1, \cdots, y_n\rp,
\textrm{sample } \tilde{y}^\lp 1\rp \sim p\lp\tilde{y}\given \theta^\lp 1\rp \rp \\
&\vdots \\
&\textrm{sample } \theta^\lp S\rp \sim p\lp\theta\given y_1, \cdots, y_n\rp,
\textrm{sample } \tilde{y}^\lp S\rp \sim p\lp\tilde{y}\given \theta^\lp S\rp \rp
\end{align}
In other words, you sample from the posterior of $\theta$ and use that sample to sample
$\tilde{Y}$. These joint samples constitute $S$ indepdendent samples from the joint posterior
distribution of $\lp \theta, \tilde{Y}\rp$, and the samples of $\tilde{Y}$ constitute
$S$ independent samples from the marginal posterior distribution of $\tilde{Y}$ which 
is the posterior predictive distribution.

# 4.4 Posterior predictive model checking

We can plot the posterior predictive distribution (made from sampling) and the 
empirical distribution of our data to see whether they match. You would imagine 
that the distributions should be similar or else it seems that the predictive
distribution would not be good for predicting future values from the empirical
distribution.

There are a couple resons the empirical (observed) distribution and predictive distribution 
may not look similar:

* the empirical distribution doesn't match the exact distribution of the population from which the data
were sampled. This could be due to sampling variability, especially when the number of samples is small
* our model is not the correct model for this data

We can evaluate these two possibilities by performing repeated simulations of the observed
dataset and checking to see how likely the observed data are given repeated samplings of the 
posterior. Note that you could check "how likely the observed data are" in different ways. 
For instance, you could compare the empirical mean versus the means from the samplings. 
In the book, they compare the ratio of 1s versus 2s for a Poisson model. The ratio for the
observed data looks much different than the ratios from the samplings which suggests the model
is not accurately capturing the shape of the distribution. However, just because a model doesn't
accurately capture one part of the data doesn't mean it's useless. The model may still recapitulate
the mean or variance. For instance, the Poisson model may still provide consistente estimates of
the mean and variance for a distribution where the mean and variance are roughly equal.

78a43860c32a3c985a756d3b4dedc3d00c2eb57a
ipynb
Jupyter Notebook
readings/A First Course in Bayesian Statistical Methods/4. Monte Carlo approximation.ipynb
rivas-lab/presentations-readings
3d9af7fba9c01ec67e80c12468106ab2cbd0ab24
MIT
2017-02-09T19:23:06.000Z
2017-03-09T18:23:37.000Z
readings/A First Course in Bayesian Statistical Methods/4. Monte Carlo approximation.ipynb
rivas-lab/reading
3d9af7fba9c01ec67e80c12468106ab2cbd0ab24
MIT
readings/A First Course in Bayesian Statistical Methods/4. Monte Carlo approximation.ipynb
rivas-lab/reading
3d9af7fba9c01ec67e80c12468106ab2cbd0ab24
MIT
Qwen/Qwen-72B
1. YES
2. YES