$\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}}$

# 3.1 The binomial model

For exchangeable binary data, our joint beliefs about the i.i.d. responses $Y_1, \cdots, Y_n$
are well approximated by

* our beliefs about $\theta = \sum_{i=1}^N Y_i / N$
* the model that, conditional on $\theta$, the $Y_i$'s are i.i.d. binary random variables 
with expectation $\theta$

If $Y_1, \cdots, Y_n \given \theta$ are i.i.d. binary$\lp\theta\rp$,

$$p\lp\theta \given y_1, \cdots, y_n \rp = \theta^{\sum y_i} \lp 1-\theta \rp^{n-\sum y_i}
\times p\lp \theta \rp / p\lp y_1, \cdots, y_n \rp$$

Given observations $y_1, \cdots, y_n$, we can use Bayes' rule to calculate 
$p\lp \theta | \given y_1, \cdots, y_n \rp$. If we have a uniform prior, the denominator $p\lp y_1, \cdots, y_n \rp$
is known from calculus to be the beta distribution. An uncertain quantity $\theta$, known 
to be between 0 and 1, has a $\textrm{beta}\lp a,b\rp$ distribution if
$$p\lp\theta\rp = \textrm{dbeta}\lp\theta, a, b\rp = 
\frac{\Gamma\lp a+b \rp}{\Gamma\lp a \rp \Gamma\lp b \rp}\theta^{a-1} \lp 1-\theta \rp^{b-1}$$

for $0 \leq \theta \leq 1$. In the case of our binary data, $a$ is the number of success (e.g. $Y_i = 1$) and 
$b$ is the number of failures (e.g. $Y_i = 0$). For a random variable with a beta distribution,

* $\textrm{mode}\l[\theta\r] = \frac{ a-1}{\lp a-1 \rp + \lp b-1 \rp}$ if $a>1$, $b>1$
* $\ev\l[\theta\r] = \frac{a}{a + b}$
* $\var\l[\theta\r] = \frac{ab}{\lp a+b+1 \rp\lp a+b\rp^2} 
= \frac{\ev\l[\theta\r] \ev\l[1 - \theta\r]}{a + b + 1}$

## 3.1.1 Inference for exchangeable binary data

If we compare the relative probabilities of any two $\theta$ values $\theta_a$ and $\theta_b$ 
by calculating $\frac{p\lp\theta_a \given y_1, \cdots, y_n\rp}{p\lp\theta_b \given y_1, \cdots, y_n\rp}$,
we will see that the resulting expression is only a function of $\sum_{i=1}^n y_i$. This means that
$\sum_{i=1}^n Y_i$ contains all the information about $\theta$ available from the data, so we say
that $\sum_{i=1}^n Y_i$ is a **sufficient statistic** for $\theta$ and $p\lp y_1, \cdots, y_n \given \theta\rp$.
Basically, the individual binary responses don't give any more info than knowing in total how many
responses were 1. It is sufficient to know $\sum_{i=1}^n Y_i$ in order to make inferences about $\theta$. 
In this case where $Y_1, \cdots, Y_n \given \theta$ are i.i.d. $\textrm{binary}\lp\theta\rp$ random variables,
the sufficient statistic $\sum_{i=1}^n Y_i$ has a binomial distribution with parameters $\lp n, \theta\rp$.

A random variable $Y\in \lb 0,1,\cdots,n\rb$ has a $\textrm{binomial}\lp n, \theta\rp$ distribution if
$$\pr\lp Y=y\given\theta\rp = \textrm{dbinom}\lp y,n,\theta\rp = {{n}\choose{y}} \theta^y\lp 1 - \theta\rp^{n-y},
y\in\lb 0,1,\cdots,n\rb$$
* $\ev\l[Y \given\theta\r] = n\theta$
* $\var\l[Y\given\theta\r] = n\theta\lp 1 - \theta\rp$

### Posterior inference under a uniform prior distribution

Having observed $Y=y$, our task is to obtain the posterior distribution of $\theta$:
$$p\lp \theta \given y\rp = \frac{p\lp y \given \theta\rp p\lp \theta \rp}{p\lp y \rp}$$
Using the binomial pdf and the fact that the denominator is the beta distribution, we 
find that 
$$p\lp\theta\given y\rp = \textrm{beta}\lp y + 1, n-y+1 \rp$$

### Posterior distributions under beta prior distributions

The uniform prior distribution has $p\lp\theta\rp = 1$ for all $\theta \in \l[0,1\r]$. 
This is equivalent to a beta distribution with $a=1,b=1$. Note that $\Gamma\lp 1 \rp = 1$
by convention. We saw above that if $\theta \sim \textrm{beta}\lp 1,1\rp \textrm{(uniform)}$
and $Y \sim \textrm{binomial}\lp n, \theta \rp$, then 
$\lb \theta \given Y=y \rp \sim \textrm{beta}\lp 1+y, 1+n-y\rp$. If we use a beta prior 
$\theta \sim \textrm{beta}\lp a,b \rp$, then
$p\lp\theta \given y\rp = \textrm{dbeta}\lp\theta, a + y, b + n - y\rp$. So the ones are 
replaced with $a$ and $b$. We can see this from
\begin{align}
p\lp \theta \given y\rp &= \frac{p\lp \theta \rp p\lp y \given \theta\rp}{p\lp y\rp} \\
&= \frac{1}{p\lp y \rp} \times \frac{\Gamma\lp a+b \rp}{\Gamma\lp a \rp \Gamma\lp b \rp}\theta^{a-1} \lp 1-\theta \rp^{b-1} \times {{n}\choose{y}} \theta^y\lp 1 - \theta\rp^{n-y}
\end{align}

The first part of this expression $\frac{1}{p\lp y \rp}$ is the normalizing factor. The second part 
comes from the fact that the prior $p\lp \theta \rp$ is a beta distribution. The third part
is the binomial pdf. We can then write
\begin{align}
p\lp \theta \given y\rp &= c\lp n, y, a, b\rp \times\theta^{a+y-1}\lp 1 - \theta\rp^{b+n-y-1} \\
&= \textrm{dbeta}\lp \theta, a + y, b + n -y\rp
\end{align}

The first line says that $p\lp \theta \given y\rp$ is proportional to 
$\theta^{a+y-1}\lp 1 - \theta\rp^{b+n-y-1}$ which has the shape of a beta distribution. 
Since $p\lp \theta \given y\rp$ and the beta density must both integrate to 1, this means
that $p\lp \theta \given y\rp$ and the beta density are in fact the same. This **trick** 
will be used throughout the book:

* We will recognize the posterior distribution is proportional to a known probability density,
and therefore must equal that density.

### Conjugacy

We have shown that a beta prior distribution and a binomial sampling model lead to a beta
posterior distribution. We say that the class of beta priors is **conjugate** for the 
binomial sampling model.

**Definition 4 (Conjugate)**: A class $\mathcal{P}$ of prior distributions for $\theta$
is called conjugate for a sampling model $p\lp y\given\theta\rp$ if 
$$p\lp \theta \rp \in \mathcal{P} \rightarrow p\lp\theta\given y\rp \in \mathcal{P}$$
In other words, the distribution is conjugate when the prior and posterior are in the same family.

While conjugate priors are easy to calculate, they may not always be the best model.
Mixtures of conjugate priors are very flexible and also computationally tractable.

### Combining information

If $\theta | Y=y \sim \textrm{beta}\lp a+y, b+n-y\rp$, then $E\l[ \theta \given y\r] = \frac{a+y}{a+b+n}$.
This **posterior expectation** is a combination of prior and data information:
\begin{align}
E\l[ \theta \given y\r] &= \frac{a+y}{a+b+n} \\
&= \frac{a+b}{a+b+n}\frac{a}{a+b} + \frac{n}{a+b+n}\frac{y}{n} \\
&=\frac{a+b}{a+b+n} \times \textrm{ prior expectation } + \frac{n}{a+b+n} \times \textrm{ data average}
\end{align}
For this model, the posterior expectation (or posterior mean) is a weighted average of the prior expectation
and the sample average, with weights proportional to $a+b$ and $n$ respectively. This leads to an
interpretation of $a$ and $b$ as "prior data":

* $a \approx$ "prior number of 1s"
* $b \approx$ "prior number of 0s"
* $a + b \approx$ "prior sample size"

If our sample size $n$ is much larger than $a+b$, then we would probably want a majority of our
information about $\theta$ to come from the data as opposed to the prior. We can see that is the
case if we look back at
$$ E\l[ \theta \given y\r] 
=\frac{a+b}{a+b+n} \times \textrm{ prior expectation } + \frac{n}{a+b+n} \times \textrm{ data average}$$
When $n >> a+b$, $E\l[ \theta \given y\r]  \approx \frac{y}{n}$, the data average.

We can also reparameterize $a$ and $b$ in terms of a "prior sample size" $n' = a + b$ and "prior frequency"
$\theta' = a / n$. This may be easier to use when we are thinking of defining our prior because we can
define what we think the frequency might be ($\theta'$) and how many samples worth of data we are sure
($n'$). In particular, we can choose $n'$ based on how many samples we know we will get so we can sort
of specify how much weight we are putting on our prior.

### Prediction

We can create a predictive distribution for data we have yet to observe.
Let's consider binary data with $y_1, \cdots, y_n$ as the observed outcomes
of $n$ binary random variables. Let $\tilde{Y} \in \lb0, 1\rb$ be an additional
outcome yet to be observed. The **predictive distributon** of $\tilde{Y}$ is
the conditional distribution of $\tilde{Y}$ given $\lb Y_1=y_1, \cdots, Y_n=y_n\rb$.
For conditionally i.i.d. binary variables this distribution can be derived from
the distribution of $\tilde{Y}$ given $\theta$ and the posterior distribution of 
$\theta$: 
\begin{align}
\pr\lp\tilde{Y} = 1 \given y_1, \cdots, y_n\rp &= \ev\l[\theta \given y_1, \cdots, y_n\r] \\
&= \frac{a + \sum_{i=1}^n y_i}{a+b+n} \\
\pr\lp\tilde{Y} = 0 \given y_1, \cdots, y_n\rp &= \frac{a + \sum_{i=1}^n \lp 1-y_i \rp}{a+b+n}
\end{align}

You can see that the predictive distribution does not rely on unknown quantities but does
rely on the observed data which makes sense.

## 3.1.2 Confidence regions

**Definition 5 (Bayesian coverage)**: An interval $\l[l\lp y\rp, u\lp y\rp\r]$, based on the
observed data $Y=y$, has 95% Bayesian coverage for $\theta$ if 
$$\pr \lp l\lp y\rp \leq \theta \leq u\lp y\rp \given Y=y\rp =0.95$$

**Definition 6 (frequentist coverage)**: A random interval $\l[l\lp Y\rp, u\lp Y\rp\r]$
has 95% frequentist coverage for $\theta$ if 
$$\pr \lp l\lp Y\rp \leq \theta \leq u\lp Y\rp \given \theta\rp =0.95$$

In a sense, the frequentist and Bayesian notions of coverage describe pre- and post-experimental
coverage respectively. 

My interpretation of this is that when you calculate a frequentist coverage region (same
as confidence interval as far as I know), the interpretation is based on *future* experiments. 
For instance, typically we think of the frequentist 95% confidence interval as having a 95%
chance of encompassing the true parameter value. However, the confidence interval actually means
that there is a 95% chance that the confidence interval from a *future* experiment will contain
the true parameter value. When we perform the experiment, the formula for the confidence interval
is specified before the experiment, so the experiment itself is the "future" experiment. The
calculation of the confidence interval is the realization of the idea that the confidence interval
calculated from this "future" experiment has a 95% chance of containing the true parameter value.

In other words, in the frequentist case, you believe that there is a "true" value of $\theta$, 
so before you even do your experiment,
you expect that the resulting confidence interval will contain the "true" value of $\theta$
with 95% probability.

The types of intervals constructed in this book are equal to the frequentist confidence intervals
asymptotically. Bayesian confidence intervals are also known as **credible regions** or credible intervals.

### Quantile-based interval

To make a $100 \times \lp 1- \alpha\rp \%$ quantile-based confidence interval,
find numbers $\theta_{\alpha / 2} < \theta_{1 - \alpha / 2}$ such that

1. $\pr\lp \theta \leq \theta_{\alpha / 2} \given Y=y\rp = \alpha / 2$
2. $\pr\lp \theta \geq \theta_{1 - \alpha / 2} \given Y=y\rp = \alpha / 2$

The numbers $\theta_{\alpha / 2}$, $\theta_{1 - \alpha / 2}$ are the $\alpha / 2$
and $1 - \alpha / 2$ posterior quantiles of $\theta$, so 
$$\pr\lp \theta \in \l[\theta_{\alpha / 2}, \theta_{1 - \alpha / 2}\r] \given Y=y \rp = 1 - \alpha$$

### Highest posterior density (HPD) region

**Definition 7 (HPD region)**: A $100 \times \lp 1-\alpha\rp\%$ HPD region consists of a 
subset of the parameter space, $s\lp y\rp \subset \Theta$ such that 

1. $\pr\lp\theta \in s\lp y\rp \given Y=y\rp = 1 - \alpha$
2. If $\theta_a \in s\lp y \rp$, and $\theta_b \notin s\lp y \rp$, then 
$p\lp \theta_a \given Y=y\rp > p\lp\theta_b \given Y=y\rp$

The first point says that the subset $s\lp y\rp$ contains $1-\alpha$ of the posterior probability.
This is the same as the $1-\alpha$ confidence interval. The second point says that every point in
$s\lp y \rp$ has a higher posterior probability than every point outside of $s\lp y \rp$. This
is not necessarily true for the $1-\alpha$ confidence interval. If the posterior is bimodal, the HPD
may not be a continuous region. You can imagine calculate the HPD by starting a horizontal line
at the top of a plot of the posterior distribution and moving it down until there is $1-\alpha$ of the
posterior above the line. This would be the HPD region.

# 3.2 The Poisson model

When our sample space is $\mathcal{Y} = \lb 0, 1, 2, \cdots \rb$, perhaps the simplest model
is the Poisson model. A random variable $Y$ has a Poisson distribution with mean $\theta$ if
$$p\lp Y=y \given \theta \rp = \textrm{dpois}\lp y,\theta\rp = \theta^ye^{-\theta}/y!$$
for $y \in \lb 0, 1, 2, \cdots \rb$. A Poisson random variable has

* $\ev \l[Y\given\theta\r] = \theta$
* $\textrm{Var}\l[Y\given\theta\r] = \theta$

The Poisson family of distributions is said to have a mean-variance relationship because
if one Poisson has a larger mean than another, it will also have a larger variance.

## 3.2.1 Posterior inference

If we model $Y_1, \cdots, Y_n$ as i.i.d. Poisson with mean $\theta$, then the joint pdf
of our sample data is
\begin{align}
\pr\lp Y_1=y_1, \cdots, Y_n=y_n \given\theta\rp &= \prod_{i=1}^n p\lp y_i\given\theta\rp \\
&=\prod_{i=1}^n\frac{1}{y_i!}\theta^{y_i}e^{-\theta} \\
&=c\lp y_1, \cdots, y_n \rp\theta^{\sum y_i}e^{-n\theta}
\end{align}
We can compare two values of $\theta$ *a posteriori* and see
\begin{align}
\frac{p\lp \theta_a \given y \rp}{p \lp \theta_b \given y \rp} &= 
\frac{e^{-n\theta_a} \theta_a^{\sum y_i} p\lp\theta_a\rp}{e^{-n\theta_b} \theta_b^{\sum y_i} p\lp\theta_b\rp}
\end{align}
As with the binary model, $\sum_{i=1}^n Y_i$ contains all the information about $\theta$ that is 
available in the data, so $\sum_{i=1}^n Y_i$ is a sufficient statistic. Furthermore,
$\lb\sum_{i=1}^n Y_i\rb \sim \textrm{Poisson}\lp n\theta\rp$.

### Conjugate prior

A class of prior densities is conjugate for the sampling model $p\lp y_1, \cdots, y_n\rp$ if the
posterior distribution is also in the class. For the Poisson sampling model, our posterior for
$\theta$ has the following form
\begin{align}
p\lp y_1, \cdots, y_n\rp &\propto p\lp \theta\rp \times p\lp y_1, \cdots, y_n \given \theta\rp
&\propto p\lp \theta\rp \times \theta^{\sum y_i} e^{-n\theta}
\end{align}
We need to find a conjugate class of densities. The simplest class of densities of this form
is known as the family of gamma distributions.

### Gamma distribution

An uncertain positive quantity $\theta$ has a $\textrm{gamma}\lp a,b \rp$ distribution if 
$$p\lp \theta \rp = \textrm{dgamma}\lp\theta,a,b\rp = \frac{b^a}{\Gamma\lp a\rp} \theta^{a-1}e^{-b\theta}$$
for $\theta,a,b > 0$. For such a random variable

* $\ev\l[\theta\r] = a/b$
* $\var\l[\theta\r] = a/b^2$
* $\textrm{mode}\l[\theta\r] = \begin{cases} \lp a-1\rp/b & \mbox{if }a>1 \\
0 & \mbox{if }a\leq 1\end{cases}$

### Posterior distribution of $\theta$

The gamma family is conjugate for the Poisson sampling model:
$$
\begin{array}{c}
\theta \sim \textrm{gamma}\lp a,b \rp \\ 
Y_1, \cdots, Y_n \given \theta \sim \textrm{Poisson}\lp\theta\rp
\end{array} 
\Rightarrow
\lb \theta \given Y_1, \cdots, Y_n\rb \sim \textrm{gamma}\lp a + \sum_{i=1}^n Y_i, b+n\rp
$$
The posterior expectation of the $\theta$ is a convex combination of the prior expectation
and the sample average
$$\ev\l[\theta\given y_1, \cdots, y_n\r] = \frac{b}{b+n}\frac{a}{b} + \frac{n}{b+n}\frac{\sum y_i}{n}$$

* $b$ is interpreted as the number of prior observations
* $a$ is interpreted as the sum of counts from $b$ prior observations

As for the binary model, the information from the data dominates for large $n$.

We can derive the posterior predictive distribution as
$$p\lp\tilde{y} | y_1, \cdots, y_n\rp =
\frac{\Gamma\lp a + \sum y_i + \tilde{y}\rp}{\Gamma\lp\tilde{y} + 1\rp\Gamma\lp a + \sum y_i\rp}
\lp\frac{b + n}{b + n + 1}\rp^{a + \sum y_i}\lp\frac{1}{b+n_1}\rp^\tilde{y}$$
for $y\in\lb 0,1,2,\cdots\rb$. This is the negative binomial distribution with parameters
$\lp a+\sum y_i, b+n\rp$ for which 

* $\ev\l[\tilde{Y} \given y_1, \cdots, y_n\r] = \frac{a + \sum y_i}{b + n} = \ev\l[\theta \given y_1, \cdots, y_n\r]$
* $\var\l[\tilde{Y} \given y_1, \cdots, y_n\r] = \frac{a + \sum y_i}{b + n}\frac{b+n+1}{b+n} =
\var\l[\theta \given y_1, \cdots, y_n\r] \times \lp b+n+1\rp = \ev\l[\tilde{Y} \given y_1, \cdots, y_n\r] \times
\frac{b+n+1}{b+n}$

The predictive variance is to some extent a measure of our posterior uncertainty about a new sample
$\tilde{Y}$ from the population. Uncertainty about $\tilde{Y}$ stems from uncertainty about the population
(uncertainty in our estimate of $\theta$) and variability in sampling. For large $n$, the uncertainty about
$\theta$ is small ($\lp b+n+1\rp/\lp b+n\rp \approx 1$) and so the predictive variability 
becomes $\theta$, the sampling variability from the Poisson. For small $n$, the uncertainty for $\tilde{Y}$
includes additional uncertainty about $\theta$.

# 3.3 Exponential families and conjugate priors

The binomial and Poisson models are both instances of one-parameter **exponential family models**.
A one-parameter exponential family model is any model whose densities can be expressed as
$$p\lp y\given \phi\rp = h\lp y\rp c\lp\phi\rp e^{\phi t\lp y\rp}$$ where $\phi$ is the unknown
parameter and $t\lp y\rp$ is the sufficient statistic.

Combining priors of the form 
$$p\lp \phi \given n_0, t_0\rp = \kappa\lp n_0, t_0\rp c\lp \phi \rp^{n_0} e^{n_0 t_0 \phi}$$
with information from $Y_1, \cdots, Y_n \sim \iid p\lp y\given \phi\rp$ yields the following
posterior distribution
\begin{align}
p\lp \phi \given y_1, \cdots, y_n\rp &\propto p\lp\phi\rp p\lp y_1, \cdots, y_n \given \phi\rp \\
&\propto c\lp \phi\rp^{n_0 + n} \exp\lb\phi\times\l[ n_0t_0 + \sum_{i=1}^n t\lp y_i \rp\r]\rb \\
&\propto p\lp \phi \given n_0 + n, n_0t_0 + n\bar{t}\lp \mathbf{y} \rp\rp
\end{align}
where $t\lp \mathbf{y} \rp = \sum t\lp y_i \rp/n$. The similarity between the posterior and prior
distributions suggests that $n_0$ can be interpreted as a "prior sample size" and $t_0$ as a 
"prior guess" of $t\lp Y\rp$. In fact
$$\ev\l[t\lp Y\rp\r] = \ev\l[\ev\l[t\lp Y\rp \given \phi\r]\r] =
\ev\l[-c'\lp \phi\rp/c\lp\phi\rp\r] = t_0$$
so $t_0$ represents the expected value of $t\lp Y\rp$. The parameter $n_0$ is a measure of how
informative the prior is. If you want to think of the relative strength of the prior versus the
data, you can think of it that the prior distribution $p\lp \phi \given n_0,t_0\rp$
contains the same information as $n_0$ independent samples from the population.

You can take the binomial and Poisson pdfs and rewrite them in the exponential family form.


```python

```

c9d71b8d6d39486f06eb7a106c1fa0dbda9393e5
ipynb
Jupyter Notebook
readings/A First Course in Bayesian Statistical Methods/3. One-parameter models.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/3. One-parameter models.ipynb
rivas-lab/reading
3d9af7fba9c01ec67e80c12468106ab2cbd0ab24
MIT
readings/A First Course in Bayesian Statistical Methods/3. One-parameter models.ipynb
rivas-lab/reading
3d9af7fba9c01ec67e80c12468106ab2cbd0ab24
MIT
Qwen/Qwen-72B
1. YES
2. YES