04-normalgamma-03-predictive.Rmd

## Predictive Distributions {#sec:NG-predictive}

```{r packages, echo=FALSE, warning=FALSE, message=FALSE,  eval=TRUE}
library(statsr)
library(ggplot2)
```


```{r tapwater, warning=FALSE, message=FALSE, echo=FALSE}
data(tapwater)
m_0 = 35; 
n_0 = 25; 
s2_0 = 156.25;
v_0 = n_0 - 1
Y = tapwater$tthm
ybar = mean(Y)
s2 = round(var(Y),1)
n = length(Y)
n_n = n_0 + n
m_n = round((n*ybar + n_0*m_0)/n_n, 1)
v_n = v_0 + n
s2_n = round( ((n-1)*s2 + v_0*s2_0 + n_0*n*(m_0 - ybar)^2/n_n)/v_n, 1)
L = qt(.025, v_n)*sqrt(s2_n/n_n) + m_n
U = qt(.975, v_n)*sqrt(s2_n/n_n) + m_n
```


In this section, we will discuss prior and posterior **predictive** distributions of the data and show how Monte Carlo sampling from the prior predictive distribution can help select hyper-parameters, while sampling from the posterior predictive distribution can be used for predicting future events or model checking.

### Prior Predictive Distribution

We can obtain the prior predictive distribution of the data from  the joint distribution of the data and the parameters  $(\mu, \sigma^2)$ or equivalently $(\mu, \phi)$, where $\phi = 1/\sigma^2$ is the precision:

**Prior:**

$$ \begin{aligned}
 \phi &\sim \textsf{Gamma}\left(\frac{v_0}{2}, \frac{v_0 s^2_0}{2} \right) \\
 \sigma^2 & = 1/\phi \\
\mu \mid \sigma^2  &\sim  \textsf{N}(m_0, \sigma^2/n_0)
\end{aligned} $$

**Sampling model:**

$$Y_i \mid \mu,\sigma^2 \iid \No(\mu, \sigma^2) $$


**Prior predictive distribution for $Y$:**

$$\begin{aligned}
p(Y) &= \iint p(Y \mid \mu,\sigma^2) p(\mu \mid \sigma^2) p(\sigma^2) d\mu \, d\sigma^2 \\
Y &\sim t(v_0, m_0, s_0^2+s_0^2/n_0)
\end{aligned}$$

By *averaging* over the possible values of the parameters from the prior distribution in the joint distribution, technically done by a double integral, we obtain the Student t as our prior predictive distribution. For those interested, details of this derivation are provided later in an optional section.
This distribution of the observables depends only on  our four  hyper-parameters from the normal-gamma family.  We can use Monte Carlo simulation  to sample from the prior predictive distribution to help elicit prior hyper-parameters as we now illustrate with the tap water example from earlier.

### Tap Water Example (continued)

A report from the city water department suggests that levels of TTHM are expected to be between 10-60 parts per billion (ppb).  Let's see how we can use this information to create an informative conjugate prior.

**Prior Mean**
First, the normal distribution and Student t distributions are symmetric around the mean or center parameter, so we will set the prior mean $\mu$ to be at the midpoint of the interval 10-60, which  would lead to  $$m_0 = (60+10)/2 = 35$$
as our prior hyper-parameter $m_0$.

**Prior Variance**
Based on the empirical rule for bell-shaped distributions, we would expect that 95% of observations are within plus or minus two standard deviations from the mean, $\pm 2\sigma$ of $\mu$.  Using this we expect that the range of the data should be approximately $4\sigma$.  Using the  values from the report, we can use this to find our prior estimate of $\sigma$, $s_0 = (60-10)/4 = 12.5$ or 
$$s_0^2 = [(60-10)/4]^2 = 156.25$$

**Prior Sample Size and Degrees of Freedom**
To complete the specification, we also need to choose the prior sample size $n_0$ and degrees of freedom $v_0$. For a sample of size $n$, the sample variance has $n-1$ degrees of freedom.  Thinking about a possible historic set of data of size $n_0$ that led to the reported interval,  we will adopt that rule to obtain the prior degrees of freedom  $v_0 = n_0 - 1$, leaving only the prior sample size to be determined. We will draw samples from the prior predictive distribution and modify $n_0$ so that the simulated data agree with our prior assumptions.

### Sampling from the Prior Predictive in `R`

The following `R` code shows a simulation from the predictive distribution with the prior sample size $n_0 = 2$. Please be careful to not confuse   the prior sample size, $n_0$, that represents the precision of our prior information with  the number of Monte Carlo simulations, $S = 10000$, that are drawn from the distributions.  These Monte Carlo samples are used to estimate quantiles of the prior predictive distribution and a large value of $S$ reduces error in the Monte Carlo approximation.

```{r predictive-TTHM}
m_0 = (60+10)/2; s2_0 = ((60-10)/4)^2;
n_0 = 2; v_0 = n_0 - 1
set.seed(1234)
S = 10000
phi = rgamma(S, v_0/2, s2_0*v_0/2)
sigma = 1/sqrt(phi)
mu = rnorm(S, mean=m_0, sd=sigma/(sqrt(n_0)))
Y = rnorm(S, mu, sigma)
quantile(Y, c(0.025,0.975))
```

Let's try to understand the code. After setting the prior hyper-parameters and random seed, we begin by simulating $\phi$ from its gamma prior distribution.  We then transform $\phi$ to calculate $\sigma$.  Using the draws of $\sigma$, we feed that into the `rnorm` function to simulate $S$ values of $\mu$ for each value of $\sigma$. The Monte Carlo draws of $\mu,\sigma$ are used to generate $S$ possible values of TTHM denoted by $Y$. In the above code we are exploiting  that all of the functions for simulating from distributions can be vectorized, i.e.  we can provide all $S$ draws of $\phi$ to the functions and get a vector result back without having to write a loop.  Finally,  we obtain the empirical quantiles from our Monte Carlo sample using the `quantile` function to approximate the actual quantiles from the prior predictive distriubtion.


This forward simulation propagates uncertainty in $\mu$ and $\sigma$ to the prior predictive distribution of the data. Calculating the sample quantiles from the samples of the prior predictive for $Y$, we see that the 95% predictive interval for TTHM includes negative values. Since TTHM cannot be negative, we can adjust $n_0$ and repeat.  Since we need a narrower interval in order to exclude zero, we can increase $n_0$ until we achieve the desired quantiles.

After some trial and error, we find that the prior sample size of 25, the empirical quantiles from the prior predictive distribution are close to the range of 10 to 60 that we were given as prior information.

```{r predictive-TTHM-best}
m_0 = (60+10)/2; s2_0 = ((60-10)/4)^2;
n_0 = 25; v_0 = n_0 - 1
set.seed(1234)
phi = rgamma(10000, v_0/2, s2_0*v_0/2)
sigma = 1/sqrt(phi)
mu = rnorm(10000, mean=m_0, sd=sigma/(sqrt(n_0)))
y = rnorm(10000, mu, sigma)
quantile(y, c(0.025,0.975))
```

Figure \@ref(fig:hist-prior) shows an estimate of the prior distribution of $\mu$ in gray and the more dispersed prior predictive distribution in TTHM in orange, obtained from the Monte Carlo samples.

```{r hist-prior, fig.align="center", fig.width=5, fig.height=3, fig.cap="Prior density", echo=FALSE}
# The palette with grey:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

nsim = length(mu)

df = data.frame(parameter = c(rep("mu", nsim), rep("Y", nsim)), x = c(mu, y))
#priorpred= ggplot(data=df, aes(x=y)) + geom_histogram(aes(x=y, y=..density..)) +
#         geom_density() + geom_density(aes(x=mu), col="blue")
ggplot(data=df, aes(x=y)) +
  geom_density(aes(x=x, colour=parameter, linetype=parameter),
               size=1.2, show.legend=FALSE) +
  stat_density(aes(x=x, colour=parameter, linetype=parameter),
               geom="line",position="identity", size=1.2) +
               xlab("TTHM (ppb)") + scale_colour_manual(values=cbPalette) +
  theme(panel.background = element_rect(fill = "transparent", colour = NA),
        legend.key = element_rect(colour = "transparent", fill = NA),
        plot.background = element_rect(fill = "transparent", colour = NA),
        legend.position = c(.75, .75),
        text = element_text(size=15))
```

Using the Monte Carlo samples, we can also estimate the prior probability of negative values of TTHM by counting the number of times the simulated values are less than zero out of the total number of simulations.

```{r negative-prior}
sum(y < 0)/length(y)  # P(Y < 0) a priori
```

With the normal prior distribution, this probability will never be zero, but may be acceptably small, so we may use the conjugate normal-gamma model for analysis.

### Posterior Predictive

We can use the same strategy to generate samples from the predictive distribution of a new measurement $Y_{n+1}$ given the observed data. In mathematical terms, the posterior predictive distribution is written as

$$Y_{n+1} \mid Y_1, \ldots, Y_n \sim \St(v_n, m_n, s^2_n (1 + 1/n_n))$$

In the code, we replace the prior hyper parameters with the posterior hyper parameters from last time.

```{r post-pred}
set.seed(1234)
phi = rgamma(10000, v_n/2, s2_n*v_n/2)
sigma = 1/sqrt(phi)
post_mu = rnorm(10000, mean=m_n, sd=sigma/(sqrt(n_n)))
pred_y =  rnorm(10000,post_mu, sigma)
quantile(pred_y, c(.025, .975))
```

Figure \@ref(fig:hist-pred) shows the Monte Carlo approximation to the prior distribution of $\mu$, and the posterior distribution of $\mu$ which is shifted to the right. The prior and posterior predictive distributions are also depicted, showing how the data have updated the prior information.

```{r hist-pred, fig.align="center", fig.width=5, fig.height=3, fig.cap="Posterior densities", echo=FALSE, message=FALSE, warning=FALSE}

nsim = length(post_mu)

df = data.frame(parameter = c(rep("prior mu", nsim),
                              rep("prior predictive Y", nsim), rep("posterior mu", nsim), rep("posterior predictive Y", nsim)), x = c(mu, y, post_mu, pred_y))

ggplot(data=df, aes(x=pred_y)) +
  geom_density(aes(x=x, colour=parameter, linetype=parameter),
               size=1.2, show.legend=FALSE) +
  stat_density(aes(x=x, colour=parameter, linetype=parameter),
     geom="line",position="identity", size=1.2) +
     xlab("TTHM (ppb)") + scale_colour_manual(values=cbPalette) +
  xlab("TTHM (ppb)") +
  scale_colour_manual(values=cbPalette) +
  theme(panel.background = element_rect(fill = "transparent", colour = NA),
        legend.key = element_rect(colour = "transparent", fill = NA),
        plot.background = element_rect(fill = "transparent", colour = NA),
        legend.position=c(.75, .75),
        text = element_text(size=15))
```

Using the Monte Carlo samples from the posterior predictive distribution, we can estimate the probability that a new TTHM sample will exceed the legal limit of 80 parts per billion, which is approximately 0.06.

```{r negative-pred}
sum(pred_y > 80)/length(pred_y)  # P(Y > 80 | data)
```

### Summary

By using Monte Carlo methods, we can obtain prior and posterior predictive distributions of the data.

* Sampling from the prior predictive distribution can help with the selection of prior hyper parameters and verify that these choices reflect the prior information that is available.

* Visualizing prior predictive distributions based on Monte Carlo simulations can help explore implications of our prior assumptions such as the choice of the hyper parameters or even assume distributions.

* If samples are incompatible with known information, such as support on positive values, we may need to modify assumptions and look at other families of prior distributions.