\chapter{Single parameter problems}
This chapter reviews some of the key results for Bayesian inference of
single parameter problems studied in MAS2903.

\section{Prior and posterior distributions}
Suppose we have data $\xvece$ which we model using the probability
(density) function $f(\xvec|\theta)$, which depends on a single
parameter~$\theta$. Once we have observed the data, $f(\xvec|\theta)$
is the \textit{likelihood function} for $\theta$ and is a function
of~$\theta$ (for fixed $\xvec$) rather than of~$\xvec$ (for
fixed~$\theta$).

Also, suppose we have prior beliefs about likely values of $\theta$
expressed by a probability (density) function $\pi(\theta)$. We can
combine both pieces of information using the following version of
Bayes Theorem. The resulting distribution for $\theta$ is called the
posterior distribution for~$\theta$ as it expresses our beliefs about
$\theta$ {\it after} seeing the data. It summarises all our current
knowledge about the parameter $\theta$.

Using Bayes Theorem, the posterior probability (density) function for
$\theta$ is
\begin{align*}
\pi(\theta|\xvec)&=\frac{\pi(\theta)\,f(\xvec|\theta)}{f(\xvec)} \\
\intertext{where}
f(\xvec)&=
\begin{cases}
\int_\Theta\pi(\theta)\,f(\xvec|\theta)\,d\theta & \text{if
$\theta$ is continuous}, \\ \\
\sum_\Theta\pi(\theta)\,f(\xvec|\theta) & \text{if
$\theta$ is discrete}.
\end{cases}
\intertext{Also, as $f(\xvec)$ is not a function of $\theta$, Bayes
  Theorem can be rewritten as}
\pi(\theta|\xvec)&\propto \pi(\theta)\times f(\xvec|\theta) \\
i.e.~~\text{posterior}&\propto\text{prior}\times\text{likelihood}.
\end{align*}

\begin{example}~\\
Table~\ref{tab:botulism} shows some data on the number of cases of
\label{ex:botulism}
foodbourne botulism in England and Wales. It is believed that cases
occur at random at a constant rate~$\theta$ in time (a Poisson
process) and so can be modelled as a random sample from a Poisson
distribution with mean~$\theta$.

\begin{table}[ht]
\begin{center}
\begin{tabular}{|lcccccccc|}
\hline
Year & 1998 & 1999 & 2000 & 2001 & 2002 & 2003 & 2004 & 2005 \\
Cases & 2 & 0 & 0 & 0 & 1 & 0 & 2 & 1\\
\hline
\end{tabular}
\caption{Number of cases of foodbourne botulism in England and Wales,
  1998--2005}
\label{tab:botulism}
\end{center}
\end{table}

An expert in the epidemiology of similar diseases gives their prior
distribution for the rate~$\theta$ as a $Ga(2,1)$ distribution, with
density
\begin{equation}
\label{eq:botulismprior}
\pi(\theta)=\theta\, e^{-\theta},\quad\theta>0, 
\end{equation}
and mean $E(\theta)=2$ and variance $Var(\theta)=2$. Determine the
posterior distribution for~$\theta$.
\end{example}

\begin{soln}
The data are observations on $X_i|\theta\sim Po(\theta)$,
$i=1,2,\ldots,8$ (independent). Therefore, the likelihood function
for~$\theta$ is
\begin{align}
\label{eq:botulismlik}
f(\xvec|\theta)&=\prod_{i=1}^8 \frac{\theta^{x_i} e^{-\theta}}{x_i!},
\quad\quad\theta>0 \notag \\
&=\frac{\theta^{2+0+\cdots+1}e^{-8\theta}}{2!\times 0!\times\cdots\times1!},\quad\quad\theta>0 \notag \\
&=\frac{\theta^6 e^{-8\theta}}{4},\quad\quad\theta>0. 
\end{align}

Bayes Theorem combines the expert opinion with the observed data, and
gives the posterior density function as

\begin{align}
\label{eq:botulismpost}
\pi(\theta|\xvec)&\propto\pi(\theta)\,f(\xvec|\theta) \notag \\
&\propto\theta e^{-\theta}\times\frac{\theta^6 e^{-8\theta}}{4},\quad\quad\theta>0 \notag  \\
&=k\,\theta^7e^{-9\theta},\quad\quad\theta>0.
\end{align}
The only continuous distribution with density of the form
$k\theta^{g-1}e^{-h\theta}$, $\theta>0$ is the $Ga(g,h)$ distribution.
Therefore, the posterior distribution must be $\theta|\xvec\sim
Ga(8,9)$.
\end{soln}

Thus the data have updated our beliefs about $\theta$ from a $Ga(2,1)$
distribution to a $Ga(8,9)$ distribution. Plots of these distributions
are given in Figure~\ref{fig:botulismpost}, and
Table~\ref{tab:botulismsum} gives a summary of the main changes
induced by incorporating the data --- a $Ga(g,h)$ distribution has
mean $g/h$, variance $g/h^2$ and mode $(g-1)/h$.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth]{botulism}
\caption{Prior (dashed) and posterior (solid) densities for $\theta$} 
\label{fig:botulismpost}
\end{center}
\end{figure}

Notice that, as the mode of the likelihood function is close to that
of the prior distribution, the information in the data is consistent
with that in the prior distribution. Also there is a reduction in
variability from the prior to the posterior distributions. The
similarity between the prior beliefs and the data has reduced the
uncertainty we have about the rate~$\theta$ at which cases occur.

\begin{table}[ht]
\begin{center}
\begin{tabular}{|l|c|c|c|}
\cline{2-4}
\multicolumn{1}{c|}{~}& Prior & Likelihood & Posterior \\
\multicolumn{1}{c|}{~}& $\eqref{eq:botulismprior}$ & $\eqref{eq:botulismlik}$ 
& $\eqref{eq:botulismpost}$ \\
\hline
$Mode(\theta)$ & 1.00 & 0.75 & 0.78 \\
$E(\theta)$    & 2.00 &  --  & 0.89 \\
$SD(\theta)$   & 1.41 &  --  & 0.31 \\
\hline
\end{tabular}
\caption{Changes in beliefs about $\theta$}
\label{tab:botulismsum}
\end{center}
\end{table}

\enlargethispage*{1cm}
\begin{example}~\\
  Consider now the general case of Example~\ref{ex:botulism}:
\label{ex:poisgam}
suppose $X_i|\theta\sim Po(\theta)$, $i=1,2,\ldots,n$ (independent)
and our prior beliefs about $\theta$ are summarised by a $Ga(g,h)$
distribution (with $g$ and $h$ known), with density
\begin{equation}
\label{eq:poprior}
\pi(\theta)=\frac{h^g\,\theta^{g-1}e^{-h\theta}}{\Gamma(g)}, 
\quad\theta>0.
\end{equation}
Determine the posterior distribution for $\theta$.

\newnotes
\begin{soln}
The likelihood function for~$\theta$ is
\begin{align}
\label{eq:polik}
f(\xvec|\theta)&=\prod_{i=1}^n \frac{\theta^{x_i} e^{-\theta}}{x_i!}, 
\quad\quad\theta>0 \notag \\
&\propto\theta^{n\bar x} e^{-n\theta}, \quad\quad\theta>0. \\
\intertext{Using Bayes Theorem, the posterior density function is}
\pi(\theta|\xvec)&\propto\pi(\theta)\,f(\xvec|\theta) \notag \\
&\propto\frac{h^g\,\theta^{g-1}e^{-h\theta}}{\Gamma(g)}\times
\theta^{n\bar x} e^{-n\theta},\quad\quad\theta>0 \notag \\
\text{i.e.}\quad\pi(\theta|\xvec)
\label{eq:popost}
&=k\theta^{g+n\bar x-1}e^{-(h+n)\theta},\quad\quad\theta>0
\end{align}
where $k$ is a constant that does not depend on $\theta$. Therefore,
the posterior density takes the form $k\theta^{G-1}e^{-H\theta}$,
$\theta>0$ and so the posterior must be a gamma distribution. Thus we have
$\theta|\xvec\sim Ga(G=g+n\bar x,H=h+n)$.
\end{soln}

%%Hack for vspace
\vspace{-2cm}
\newpage 

Summary: 

If we have a random sample from a $Po(\theta)$ distribution and our
prior beliefs about $\theta$ follow a $Ga(g,h)$ distribution then,
after incorporating the data, our (posterior) beliefs about $\theta$
follow a $Ga(g+n\bar x,h+n)$ distribution.

The changes in our beliefs about $\theta$ are summarised in
Table~\ref{tab:gam}, taking $g\geq 1$.
\begin{table}[ht]
\bigskip
\begin{center}
\begin{tabular}{|l|c|c|c|}
\cline{2-4}
\multicolumn{1}{c|}{~}& Prior & Likelihood & Posterior \\
\multicolumn{1}{c|}{~}& $\eqref{eq:poprior}$ & $\eqref{eq:polik}$ & $\eqref{eq:popost}$ \\
\hline
$Mode(\theta)$ & $(g-1)/h$ & $\bar x$ & $(g+n\bar x-1)/(h+n)$ \\
$E(\theta)$ & $g/h$ & -- & $(g+n\bar x)/(h+n)$ \\
$SD(\theta)$ & $\sqrt{g}/h$ & -- & $\sqrt{g+n\bar x}/(h+n)$ \\
\hline
\end{tabular}
\caption{Changes in beliefs about $\theta$}
\label{tab:gam}
\end{center}
\end{table}
Notice that the posterior mean is greater than the prior mean if and
only if the likelihood mode is greater than the prior mean, that is,
\begin{equation*}
E(\theta|\xvec)>E(\theta)\quad\iff\quad 
Mode_\theta\{f(\xvec|\theta)\}>E(\theta).
\end{equation*}
The standard deviation of the posterior distribution is smaller than
that of the prior distribution if and only if the sample mean is not
too large, that is
\[
SD(\theta|\xvec)<SD(\theta)\quad\iff\quad 
Mode_\theta\{f(\xvec|\theta)\} < \left(2+\frac{n}{h}\right)E(\theta),
\]
and this will be true in large samples.
\end{example}

\begin{example}
  \label{ex:normal}~\\
  Suppose we have a random sample from a normal distribution. In
  Bayesian statistics, when dealing with the normal distribution, the
  mathematics is more straightforward working with the precision
  ($=1/\text{variance}$) of the distribution rather than the
  variance itself. So we will assume that this population has unknown
  mean~$\mu$ but known precision~$\tau$: $X_i|\mu\sim N(\mu,1/\tau)$,
  $i=1,2,\ldots,n$ (independent), where $\tau$ is known. Suppose our
  prior beliefs about $\mu$ can be summarised by a $N(b,1/d)$
  distribution, with probability density function
\begin{equation}
\label{eq:p11}
\pi(\mu)=
\left(\frac{d}{2\pi}\right)^{1/2}\exp\left\{-\frac{d}{2}(\mu-b)^2\right\}.
\end{equation}
Determine the posterior distribution for $\mu$. 

Hint:
\[
d(\mu-b)^2+n\tau(\bar x-\mu)^2
=(d+n\tau)\left\{\mu
-\left(\frac{db+n\tau\bar x}{d+n\tau}\right)\right\}^2+c
\]
where $c$ does not depend on $\mu$.

\begin{soln}
The likelihood function for~$\mu$ is
\begin{align*}
f(\xvec|\mu)&=\prod_{i=1}^n 
\left(\frac{\tau}{2\pi}\right)^{1/2}
\exp\left\{-\frac{\tau}{2}(x_i-\mu)^2\right\} \\
&= \left(\frac{\tau}{2\pi}\right)^{n/2}
\exp\left\{-\frac{\tau}{2}\sum_{i=1}^n (x_i-\mu)^2\right\}.
\intertext{Now}
\sum_{i=1}^n (x_i-\mu)^2
&=\sum_{i=1}^n (x_i-\bar x+\bar x-\mu)^2\\
&=\sum_{i=1}^n \{(x_i-\bar x)^2+(\bar x-\mu)^2+2(x_i-\bar x)(\bar x-\mu)\}\\
&=\sum_{i=1}^n \{(x_i-\bar x)^2+(\bar x-\mu)^2\}
+2(\bar x-\mu)\sum_{i=1}^n (x_i-\bar x)\\
&=\sum_{i=1}^n (x_i-\bar x)^2+n(\bar x-\mu)^2.
\intertext{Let $s^2=\dfrac{1}{n}\displaystyle\sum_{i=1}^n (x_i-\bar
  x)^2$ and so}
\sum_{i=1}^n (x_i-\mu)^2&=n\left[s^2+(\bar x-\mu)^2\right].
\end{align*}
Therefore
\begin{align}
\label{eq:p10}
f(\xvec|\mu)&= \left(\frac{\tau}{2\pi}\right)^{n/2}
\exp\left\{-\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\}.
\end{align}
Using Bayes Theorem, the posterior density function is, for
$\mu\in\R$
\begin{align*}
\pi(\mu|\xvec)&\propto\pi(\mu)\,f(\xvec|\mu)\\
&\propto \left(\frac{d}{2\pi}\right)^{1/2}\exp\left\{-\frac{d}{2}(\mu-b)^2\right\} \\
&\hspace*{3cm}\times \left(\frac{\tau}{2\pi}\right)^{n/2}
\exp\left\{-\frac{n\tau}{2}\left[s^2+(\bar x-\mu)^2\right]\right\} \\
&\propto\exp\left\{-\frac{1}{2}\left[d(\mu-b)^2+n\tau(\bar x-\mu)^2\right]\right\}\\
&\propto\exp\left\{-\frac{1}{2}\left[(d+n\tau)\left\{\mu
-\left(\frac{db+n\tau\bar x}{d+n\tau}\right)\right\}^2+c\right]\right\}\\
\intertext{\hfill using the hint}
&\propto\exp\left\{-\frac{1}{2}\left[(d+n\tau)\left\{\mu
-\left(\frac{db+n\tau\bar x}{d+n\tau}\right)\right\}^2\right]\right\}
\end{align*}
as $c$ does not depend on $\mu$. Let
\begin{equation}
\label{eq:p12}
B=\frac{db+n\tau\bar x}{d+n\tau}\quad\quad\text{and}\quad\quad
D=d+n\tau.
\end{equation}
Then
\begin{equation}
\label{eq:p13}
\pi(\mu|\xvec)=k\exp\left\{-\frac{D}{2}(\mu-B)^2\right\},
\end{equation}
where $k$ is a constant that does not depend on $\mu$. Therefore, the
posterior density takes the form $k\exp\{-D(\mu-B)^2/2\}$, $\mu\in\R$
and so the posterior distribution must be a normal distribution: we
have $\mu|\xvec\sim N(B,1/D)$.
\end{soln}

%%Hack for vspace
\vspace{-3cm}
\newpage
Summary: 

If we have a random sample from a $N(\mu,1/\tau)$ distribution (with
$\tau$ known) and our prior beliefs about $\mu$ follow a $N(b,1/d)$
distribution then, after incorporating the data, our (posterior)
beliefs about $\mu$ follow a $N(B,1/D)$ distribution.

The changes in our beliefs about $\mu$ are summarised in
Table~\ref{tab:norknown}. 
\begin{table}[ht]
\bigskip
\begin{center}
\begin{tabular}{|l|c|c|c|}
\cline{2-4}
\multicolumn{1}{c|}{~}& Prior & Likelihood & Posterior \\
\multicolumn{1}{c|}{~}& $\eqref{eq:p11}$ & $\eqref{eq:p10}$ & $\eqref{eq:p13}$ \\
\hline
$Mode(\mu)$ & $b$ & $\bar x$ & $(db+n\tau\bar x)/(d+n\tau)$ \\
$E(\mu)$ & $b$ & -- & $(db+n\tau\bar x)/(d+n\tau)$ \\
$Precision(\mu)$ & $d$ & -- & $d+n\tau$ \\
\hline
\end{tabular}
\caption{Changes in beliefs about $\mu$}
\label{tab:norknown}
\end{center}
\end{table}
Notice that the posterior mean is greater than the prior mean if and
only if the likelihood mode (sample mean) is greater than the prior
mean, that is
\begin{equation*}
E(\mu|\xvec)>E(\mu)\quad\iff\quad Mode_\mu\{f(\xvec|\mu)\}>E(\mu).
\end{equation*}
Also, the standard deviation of the posterior distribution is smaller
than that of the prior distribution.
\end{example}

\begin{example}~\\
The 18th century physicist Henry Cavendish made 23 experimental
\label{ex:cavendish}
  determinations of the earth's density, and these data (in $g/cm^3$)
  are given below.
\begin{center}
\begin{tabular}{|cccccccc|}
\hline
5.36&5.29&5.58&5.65&5.57&5.53&5.62&5.29\\
5.44&5.34&5.79&5.10&5.27&5.39&5.42&5.47\\
5.63&5.34&5.46&5.30&5.78&5.68&5.85&\\
\hline
\end{tabular}
\end{center}
Suppose that Cavendish asserts that the error standard deviation of
these measurements is $0.2\,g/cm^3$, and assume that they are normally
distributed with mean equal to the true earth density~$\mu$.  Using a
normal prior distribution for $\mu$ with mean $5.41\,g/cm^3$ and
standard deviation $0.4\,g/cm^3$, derive the posterior distribution
for $\mu$.

%\newstud
\begin{soln}
  From the data we calculate $\bar x=5.4848$ and $s=0.1882$.
  Therefore, the assumed standard deviation $\sigma=0.2$ is probably
  okay. We also have $\tau=1/0.2^2$, $b=5.41$, $d=1/0.4^2$ and $n=23$.
  Therefore, using Example~\ref{ex:normal}, the posterior distribution
  is $\mu|\xvec\sim N(B,1/D)$, where
\begin{align*}
B&=\frac{db+n\tau\bar x}{d+n\tau}
=\frac{5.41/0.4^2+23\times 5.4848/0.2^2}{1/0.4^2+23/0.2^2}
=5.4840
\intertext{and}
D&=d+n\tau=\frac{1}{0.4^2}+\frac{23}{0.2^2}=\frac{1}{0.0415^2}.
\end{align*}
Therefore the posterior distribution is $\mu|\xvec\sim
N(5.484,0.0415^2)$ and is shown in Figure~\ref{fig:cavendishplot}.
\end{soln}

%%Hack for vspace
\vspace{-3.5cm}

The actual mean density of the earth is $5.515\,g/cm^3$ (Wikipedia).
We can determine the (posterior) probability that the mean density is
within 0.1 of this value as follows. The posterior distribution is
$\mu|\xvec\sim N(5.484,0.0415^2)$ and so
\[
Pr(5.415<\mu<5.615|\xvec)=0.9510,
\]
calculated using the \verb$R$ command
\verb$pnorm(5.615,5.484,0.0415)-pnorm(5.415,5.484,0.0415)$.

Without the data, the only basis for determining the earth's density
is via the prior distribution. Here the prior distribution is $\mu\sim
N(5.4,0.4^2)$ and so the (prior) probability that the mean density is
within 0.2 of the (now known) true value is
\[
Pr(5.315<\mu<5.715)=0.1896,
\]
calculated using the \verb$R$ command
\verb$pnorm(5.615,5.4,0.4)-pnorm(5.415,5.4,0.4)$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth]{cavendish}
\caption{Prior (dashed) and posterior (solid) densities for the
  earth's density}
\label{fig:cavendishplot}
\end{center}
\end{figure}
\end{example}

\section{Different levels of prior knowledge}
\subsection{Substantial prior knowledge}
We have {\it substantial prior information} for $\theta$ when the
prior distribution {\it dominates} the posterior distribution, that is
$\pi(\theta|\xvec)\sim\pi(\theta)$. 

When we have substantial prior information there can be some
difficulties:
\begin{enumerate}
\item the intractability of the mathematics in deriving the posterior
distribution --- though with modern computing facilities this is less
of a problem,
\item the practical formulation of the prior distribution ---
coherently specifying prior beliefs in the form of a probability
distribution is far from straightforward.
\end{enumerate}

\subsection{Limited prior knowledge}
When prior information about $\theta$ is limited, the pragmatic
approach is to choose a distribution which makes the Bayes updating
from prior to posterior mathematically straightforward, and use what
prior information is available to determine the parameters of this
distribution. For example
\begin{itemize}
\item Poisson random sample, Gamma prior distribution
  $\longrightarrow$ Gamma posterior distribution
\item Normal random sample (known variance), Normal prior distribution
  $\longrightarrow$ Normal posterior distribution
\end{itemize}

In these examples, the prior distribution and the posterior
distribution come from the same family. This leads us to the following
definition.

\begin{definition}~\\
Suppose that data $\xvec$ are to be observed with distribution
$f(\xvec|\theta)$. A family $\mathfrak{F}$ of prior distributions for
$\theta$ is said to be {\it conjugate} to $f(\xvec|\theta)$ if
for every prior distribution $\pi(\theta)\in\mathfrak{F}$, the
posterior distribution $\pi(\theta|\xvec)$ is also in $\mathfrak{F}$.
\end{definition}

Notice that the conjugate family depends crucially on the model chosen
for the data $\xvec$. For example, the only family conjugate to the
model ``random sample from a Poisson distribution'' is the Gamma
family.

\subsection{Vague prior knowledge}
If we have very little or no prior information about the model
parameters $\theta$, we must still choose a prior distribution in
order to operate Bayes Theorem. Obviously, it would be sensible to
choose a prior distribution which is not concentrated about any
particular value, that is, one with a very large variance. In
particular, most of the information about $\theta$ will be passed
through to the posterior distribution via the data, and so we have
$\pi(\theta|\xvec)\sim f(\xvec|\theta)$.

We represent vague prior knowledge by using a prior distribution which
is conjugate to the model for $\xvec$ and which is as diffuse as
possible, that is, has as large a variance as possible.

\newnotes
\begin{example}~\\
Suppose we have a random sample from a $N(\mu,1/\tau)$ distribution
(with $\tau$ known). Determine the posterior distribution assuming a
vague prior for $\mu$.

\begin{soln}
The conjugate prior distribution is a normal distribution. We have
already seen that if the prior is $\mu\sim N(b,1/d)$ then the
posterior distribution is $\mu|\xvec\sim N(B,1/D)$ where
\begin{equation*}
B=\frac{db+n\tau\bar x}{d+n\tau}\qquad\text{and}\qquad
D=d+n\tau.
\end{equation*}
If we now make our prior knowledge vague about $\mu$ by letting the
prior variance tend to infinity ($d\to 0$), we obtain
\[
B\to\bar x\qquad\text{and}\qquad D\to n\tau.
\]
Therefore, assuming vague prior knowledge for $\mu$ results in a
$N\{\bar x,1/(n\tau)\}$ posterior distribution. 

Notice that the posterior mean is the sample mean (the likelihood
mode) and that the posterior variance $1/(n\tau)\to 0$ as
$n\to\infty$.
\end{soln}
\end{example}

%%Hack for vspace
\newpage

\begin{example}~\\
Suppose we have a random sample from a Poisson distribution, that
\label{ex:vaguegamma}
is, $X_i|\theta\sim Po(\theta)$, $i=1,2,\ldots,n$
(independent). Determine the posterior distribution assuming a vague
prior for $\theta$.
\end{example}

\textbf{\large Solution}\\~\\\noindent
  The conjugate prior distribution is a Gamma distribution. Recall
  that a $Ga(g,h)$ distribution has mean $m=g/h$ and variance
  $v=g/h^2$. Rearranging these formulae we obtain
\begin{equation*}
g=\frac{m^2}{v}\quad\quad\text{and}\quad\quad h=\frac{m}{v}. 
\end{equation*}
Clearly $g\to 0$ and $h\to 0$ as $v\to\infty$ (for fixed~$m$).  We
have seen how taking a $Ga(g,h)$ prior distribution results in a
$Ga(g+n\bar x,h+n)$ posterior distribution. Therefore, taking a vague
prior distribution will give a $Ga(n\bar x,n)$ posterior distribution.

Note that the posterior mean is $\bar x$ (the likelihood mode) and
that the posterior variance $\bar x/n\to 0$ and $n\to\infty$.

\section{Asymptotic posterior distribution}
If we have a statistical model $f(\xvec|\theta)$ for data $\xvece$,
together with a prior distribution $\pi(\theta)$ for $\theta$ then
\begin{align*}
  &\sqrt{J(\thetahat)}~(\theta-\thetahat)|\xvec \stackrel{\cal D}
  \longrightarrow N(0,1)
  \quad\quad\text{as }n\to\infty, \\
  \intertext{where $\thetahat$ is the likelihood mode and $J(\theta)$
    is the observed information}
  &J(\theta)=-\frac{\partial^2}{\partial\theta^2}\,
  \log f(\xvec|\theta).
\end{align*}

This means that, with increasing amounts of data, the posterior
distribution looks more and more like a normal distribution. The
result also gives us a useful approximation to the
posterior distribution for $\theta$ when $n$ is large:
\[
\theta|\xvec\sim N\{\hat\theta,J(\hat\theta)^{-1}\}
\quad\quad\text{approximately}.
\]
Note that this limiting result is similar to one used in Frequentist
statistics for the distribution of the maximum likelihood estimator,
namely
\begin{align*}
  &\sqrt{I(\theta)}~(\thetahat-\theta)\stackrel{\cal D}
  \longrightarrow N(0,1)
  \quad\quad\text{as }n\to\infty,
\end{align*}
where Fisher's information $I(\theta)$ is the expected value of the
observed information, where the expectation is taken over the
distribution of $\Xvec|\theta$, that is,
$I(\theta)=E_{\Xvec|\theta}[J(\theta)]$. You may also have seen this
result written as an approximation to the distribution of the maximum
likelihood estimator in large samples, namely
\[
\thetahat\sim N\{\theta,I(\theta)^{-1}\}
\quad\quad\text{approximately}.
\]

\begin{example}~\\
Suppose we have a random sample from a $N(\mu,1/\tau)$ distribution
(with $\tau$ known). Determine the asymptotic posterior distribution
for~$\mu$. 

Recall that
\[
f(\xvec|\mu)=
\left(\frac{\tau}{2\pi}\right)^{n/2}
\exp\left\{-\frac{\tau}{2}\sum_{i=1}^n (x_i-\mu)^2\right\},
\]
and therefore 
\begin{align*}
\log\,&f(\xvec|\mu)=\frac{n}{2}\log\tau-\frac{n}{2}\log(2\pi)
-\frac{\tau}{2}\sum_{i=1}^n (x_i-\mu)^2 \\
&\Rightarrow\quad \frac{\partial}{\partial\mu} \log
f(\xvec|\mu)=-\frac{\tau}{2}\times\sum_{i=1}^n -2(x_i-\mu)
=\tau\sum_{i=1}^n (x_i-\mu) =n\tau(\bar x-\mu)\\
&\Rightarrow\quad \frac{\partial^2}{\partial\mu^2} \log
f(\xvec|\mu)=-n\tau
\qquad\Rightarrow\quad J(\mu)=-\frac{\partial^2}{\partial\mu^2} \log
f(\xvec|\mu)=n\tau.
\end{align*}


\begin{soln}
We have
\begin{align*}
\frac{\partial}{\partial\mu} \log f(\xvec|\mu)=0
\quad\quad&\Longrightarrow\quad\quad
\hat\mu=\bar x \\
&\Longrightarrow\quad\quad
J(\hat\mu)=n\tau \\
&\Longrightarrow\quad\quad
J(\hat\mu)^{-1}=\frac{1}{n\tau}.
\end{align*}
Therefore, for large $n$, the (approximate) posterior distribution for
$\mu$ is
\[
\mu|\xvec\sim N\left(\bar x,\,\frac{1}{n\tau}\right).
\]
\end{soln}
Here the asymptotic posterior distribution is the same as the
posterior distribution under vague prior knowledge.
\end{example}


\section{Bayesian inference}

The posterior distribution $\pi(\thetavec|\xvec)$ summarises all our
information about $\thetavec$ to date. However, sometimes it is
helpful to reduce this distribution to a few key summary measures.

\subsection{Estimation}
\label{sec:univariate-estimation}
\subsubsection*{Point estimates}
There are many useful summaries for a typical value of a random
variable with a particular distribution; for example, the mean, mode
and median. The mode is used more often as a summary than is the case
in frequentist statistics.

\subsubsection*{Confidence intervals/regions}
A more useful summary of the posterior distribution is one which also
reflects its variation. For example, a $100(1-\alpha)\%$ {\it Bayesian
  confidence interval} for $\theta$ is any region $C_\alpha$ that
satisfies $Pr(\theta\in C_\alpha|\xvec)=1-\alpha$.  If $\theta$
is a continuous quantity with posterior probability density function
$\pi(\theta|\xvec)$ then
\begin{equation*}
\int_{C_\alpha} \pi(\theta|\xvec)\,d\theta = 1-\alpha.
\end{equation*}
The usual correction is made for discrete $\theta$, that is, we
take the largest region $C_\alpha$ such that $Pr(\theta\in
C_\alpha|\xvec)\leq 1-\alpha$. Bayesian confidence intervals are
sometimes called {\it credible regions} or {\it plausible regions}.
Clearly these intervals are not unique, since there will be many
intervals with the correct probability coverage for a given posterior
distribution.

A $100(1-\alpha)\%$ \textit{highest density interval} (HDI) for
$\theta$ is the region
\[
C_\alpha=\{\theta:~\pi(\theta|\xvec)\geq\gamma\}
\] 
where $\gamma$ is chosen so that $Pr(\theta\in
C_\alpha|\xvec)=1-\alpha$.  This region is sometimes called a
\textit{most plausible Bayesian confidence interval}. If the posterior
distribution has many modes then it is possible that the HDI will be
the union of several disjoint regions.  Also, if the posterior
distribution is unimodal (has one mode) and symmetric about its mean
then the HDI is an equi-tailed interval, that is, takes the form
$C_\alpha=(a,b)$, where
$Pr(\theta<a|\xvec)=Pr(\theta>b|\xvec)=\alpha/2$; see Figure~\ref{fig:hdi}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth]{hdi}
\caption{Construction of an HDI for a symmetric posterior density}
\label{fig:hdi}
\end{center}
\end{figure}

%%Hack for vspace
\newpage
\subsubsection*{Interpretation of confidence intervals/regions}
Suppose $C_B$ is a 95\% Bayesian confidence interval for $\theta$
and $C_F$ is a 95\% frequentist confidence interval for $\theta$.
These intervals do not have the same interpretation:
\begin{myitemize}
\item the probability that $C_B$ contains $\theta$ is 0.95;
\item the probability that $C_F$ contains $\theta$ is either 0 or 1
--- since $\theta$ does not have a (non-degenerate) probability
distribution;
\item the interval $C_F$ covers the true value $\theta$ on 95\% of
  occasions --- in repeated applications of the formula.
\end{myitemize}

\begin{example}~\\
Suppose we have a random sample $\xvece$ from a $N(\mu,1/\tau)$
distribution (where $\tau$ is known). We have seen that, assuming
vague prior knowledge, the posterior distribution is $\mu|\xvec\sim
N\{\bar x,1/(n\tau)\}$. Determine the $100(1-\alpha)\%$ HDI for~$\mu$.

\begin{soln}
  This distribution has a symmetric bell shape and so the HDI is an
  equi-tailed interval $C_\alpha=(a,b)$ with
  $Pr(\mu<a|\xvec)=\alpha/2$ and $Pr(\mu>b|\xvec)=\alpha/2$, that is,
\[
a=\bar x - \frac{z_{\alpha/2}}{\sqrt{n\tau}}
\quad\quad\text{and}\quad\quad
b=\bar x + \frac{z_{\alpha/2}}{\sqrt{n\tau}},
\]
where $z_\alpha$ is the upper $\alpha$-quantile of the $N(0,1)$
distribution. For example, the 95\% HDI for~$\mu$ is 
\[
\left(\bar x - \frac{1.96}{\sqrt{n\tau}},\,
\bar x + \frac{1.96}{\sqrt{n\tau}}\right).
\]
\end{soln}
Note that this interval is numerically identical to the 95\%
frequentist confidence interval for the (population) mean of a normal
random sample with known variance. However, the interpretation is very
different.
\end{example}

\begin{example}~\\
  Recall Example~\ref{ex:botulism} on the number of cases of
  foodbourne botulism in England and Wales. The data were modelled as
  a random sample from a Poisson distribution with mean~$\theta$.
  Using a $Ga(2,1)$ prior distribution, we found the posterior
  distribution to be $\theta|\xvec\sim Ga(8,9)$. This posterior
  density is shown in Figure~\ref{fig:botulismpost2}. Determine the 
  $100(1-\alpha)\%$ HDI for~$\theta$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth]{botulism2}
\caption{Posterior density for $\theta$} 
\label{fig:botulismpost2}
\end{center}
\end{figure}

\begin{soln}
  The HDI must take the form $C_\alpha=(a,b)$ if it is to include the
  values of $\theta$ with the highest probability density. Suppose
  that $F(\cdot)$ and $f(\cdot)$ are the posterior distribution and density
  functions. Then the end-points $a$ and $b$ must satisfy
\begin{align*}
Pr(a<\theta<b|\xvec)&=F(b)-F(a)=1-\alpha\\
\intertext{and}
\pi(\theta=a|\xvec)-\pi(\theta=b|\xvec)&=f(a)-f(b)=0.
\end{align*}
Unfortunately, there is no simple analytical solution to these
equations and so numerical methods have to be employed to determine
$a$ and $b$. However, if we have the quantile function $F^{-1}(\cdot)$
for $\theta|\xvec$ then we can find the solution by noticing that we
can write $b$ in terms of $a$:
\[
b=F^{-1}\{1-\alpha+F(a)\},
\]
for $0<a<F^{-1}(\alpha)$. Therefore, we can determine the correct choice
of $a$ by minimising the function
\[
g(a)=\bigl(f(a)-f[F^{-1}\{1-\alpha+F(a)\}]\bigr)^2.
\]
within the range $0<a<F^{-1}(\alpha)$. The values of $a$ and $b$ can
be determined using the optimizer function \texttt{optimize} in
\texttt{R}.
\end{soln}

%%Hack for vspace
\newpage

The \texttt{R} package \texttt{nclbayes} contains functions to
determine the HDI for several distributions. The function for the
Gamma distribution is \texttt{hdiGamma} and we can calculate the 95\%
HDI for the $Ga(8,9)$ posterior distribution by using the commands
\begin{verbatim}
library(nclbayes)
hdiGamma(p=0.95,a=8,b=9)
\end{verbatim}

Taking $1-\alpha=0.95$ and using such \texttt{R} code gives
$a=0.3304362$ and $b=1.5146208$. To check this answer, \texttt{R}
gives $Pr(a<\theta<b|\xvec)=0.95$, $\pi(\theta=b|\xvec)=0.1877215$ and
$\pi(\theta=a|\xvec)=0.1877427$.  Thus the 95\% HDI is
$(0.3304362,1.514621)$.


The package also has functions \texttt{hdiBeta} for the Beta
distribution and \texttt{hdiInvchi} for the Inv-Chi distribution
(introduced in Chapter~2).
\end{example}

\subsection{Prediction}
Much of statistical inference (both Frequentist and Bayesian) is aimed
towards making statements about a parameter $\theta$. Often the
inferences are used as a yardstick for similar future experiments. For
example, we may want to predict the outcome when the experiment is
performed again.

Clearly there will be uncertainty about the future outcome of an
experiment. Suppose this future outcome $Y$ is described by a
probability (density) function $f(y|\theta)$. There are several
ways we could make inferences about what values of $Y$ are likely. For
example, if we have an estimate $\thetahat$ of $\theta$ we might
base our inferences on $f(y|\theta=\thetahat)$. Obviously this
is not the best we can do, as such inferences ignore the fact that it
is very unlikely that $\theta=\thetahat$. 

Implicit in the Bayesian framework is the concept of the
\textit{predictive distribution}. This distribution describes how
likely are different outcomes of a future experiment. The predictive
probability (density) function is calculated as
\[
f(y|\xvec)=\int_\Theta
f(y|\theta)\,\pi(\theta|\xvec)\,d\theta
\]
when $\theta$ is a continuous quantity. From this equation, we can
see that the predictive distribution is formed by weighting the
possible values of $\theta$ in the future experiment
$f(y|\theta)$ by how likely we believe they are to occur
$\pi(\theta|\xvec)$.

If the true value of $\theta$ were known, say $\theta_0$, then
any prediction can do no better than one based on
$f(y|\theta=\theta_0)$. However, as (generally) $\theta$ is
unknown, the predictive distribution is used as the next best
alternative.

We can use the predictive distribution to provide a useful range of
plausible values for the outcome of a future experiment. This {\it
prediction interval} is similar to a HDI interval.  A
$100(1-\alpha)\%$ {\it prediction interval} for $Y$ is the region
$C_\alpha=\{y:~f(y|\xvec)\geq\gamma\}$ where $\gamma$ is chosen so
that $Pr(Y\in C_\alpha|\xvec)=1-\alpha$.

\begin{example}~\\
  Recall Example~\ref{ex:botulism} on the number of cases of
\label{ex:botulism2}
foodbourne botulism in England and Wales. The data for 1998--2005 were
modelled by a Poisson distribution with mean~$\theta$. Using a
$Ga(2,1)$ prior distribution, we found the posterior distribution to
be $\theta|\xvec\sim Ga(8,9)$. Determine the predictive distribution
for the number of cases for the following year (2006).

\begin{soln}
  Suppose the number of cases in 2006 is $Y$, with $Y|\theta\sim
  Po(\theta)$. The predictive probability function of $Y$ is, for
  $y=0,1,\ldots$
\begin{align*}
f(y|\xvec)
&=\int_\Theta f(y|\theta)\,\pi(\theta|\xvec)\,d\theta \\
&=\int_0^\infty \frac{\theta^ye^{-\theta}}{y!}
\times \frac{9^8\theta^7e^{-9\theta}}{\Gamma(8)}\,d\theta
\end{align*}
\begin{align*}
\phantom{f(y|\xvec)}
&=\frac{9^8}{y!\Gamma(8)}
\int_0^\infty \theta^{y+7}e^{-10\theta}\,d\theta \\
&=\frac{9^8}{y!\Gamma(8)}\times\frac{\Gamma(y+8)}{10^{y+8}}\\
&=\frac{(y+7)!}{y!7!}\times 0.9^8\times 0.1^y\\
&=\binom{y+7}{7}\times 0.9^8\times 0.1^y.
\end{align*}
\end{soln}
You may not recognise this probability function but it is related to
that of a negative binomial distribution. Suppose $Z\sim NegBin(r,p)$
with probability function
\[
Pr(Z=z)=\binom{z-1}{r-1}p^r(1-p)^{z-r},\quad z=r,r+1,\ldots.
\]
Then $W=Z-r$ has probability function
\[
Pr(W=w)=Pr(Z=w+r)=\binom{w+r-1}{r-1}p^r(1-p)^w,\quad w=0,1,\ldots.
\]
This is the same probability function as our predictive probability
function, with $r=8$ and $p=0.9$.  Therefore $Y|\xvec\sim
NegBin(8,0.9)-8$.  Note that, unfortunately \texttt{R} also calls the
distribution of~$W$ a negative binomial distribution with
parameters~$r$ and~$p$. To distinguish between this distribution and
the $NegBin(r,p)$ distribution used above, we shall denote the
distribution of~$W$ as a $NegBin_\texttt{R}(r,p)$ distribution -- it
has mean $r(1-p)/p$ and variance $r(1-p)/p^2$. Thus $Y|\xvec\sim
NegBin_\texttt{R}(8,0.9)$.
\end{example}

We can compare this predictive distribution with a naive predictive
distribution based on an estimate of $\theta$. Here we shall base our
naive predictive distribution on the maximum likelihood estimate
$\thetahat=0.75$, that is, use the distribution
$Y|\theta=\hat\theta\sim Po(0.75)$. Thus, the naive predictive
probability function is
\[
f(y|\theta=\thetahat)=\frac{0.75^y\,e^{-0.75}}{y!},\quad y=0,1,\ldots~.
\]

Numerical values for the predictive and naive predictive probability
functions are given in Table~\ref{tab:predbetabin}.
\begin{table}[ht]
  \bigskip
\begin{center}
\begin{tabular}{|r|c|c|}
\cline{2-3}
\multicolumn{1}{c|}{~}& correct & naive \\
\cline{1-1}
$y$ & $f(y|\xvec)$ & $f(y|\theta=\thetahat)$ \\
\hline
0 & 0.430 & 0.472 \\
1 & 0.344 & 0.354 \\
2 & 0.155 & 0.133 \\
3 & 0.052 & 0.033 \\
4 & 0.014 & 0.006 \\
5 & 0.003 & 0.001 \\
$\geq 6$ & 0.005 & 0.002 \\
\hline
\end{tabular}
\end{center}
\caption{Predictive and naive predictive probability functions}
\label{tab:predbetabin}
\end{table}

Again, the naive predictive distribution is a predictive distribution
which, instead of using the correct posterior distribution, uses a
degenerate posterior distribution $\pi^*(\theta|\xvec)$ which
essentially allows only one value: $Pr_{\pi^*}(\theta=0.75|\xvec)=1$
and standard deviation $SD_{\pi^*}(\theta|\xvec)=0$. Note that the
correct posterior standard deviation of $\theta$ is
$SD_\pi(\theta|\xvec)=\sqrt{8}/9=0.314$. Using a degenerate posterior
distribution results in the naive predictive distribution having too
small a standard deviation:
\begin{equation*}
SD(Y|x=1)=
\begin{cases} 0.994 & \text{using the correct }\pi(\theta|\xvec) \\
              0.866 & \text{using the naive }\pi^*(\theta|\xvec), \\
\end{cases}
\end{equation*}
these values being calculated from 
$NegBin_\texttt{R}(8,0.9)$ and $Po(0.75)$ distributions.

Using the numerical table of predictive probabilities, we can see that
$\{0,1,2\}$ is a 92.9\% prediction set/interval. This is to be
contrasted with the more ``optimistic'' calculation using the naive
predictive distribution which shows that $\{0,1,2\}$ is a 95.9\%
prediction set/interval.

\subsection*{Candidate's formula}
In the previous example, a non-trivial integral had to be evaluated.
However, when the past data $\xvec$ and future data $y$ are
independent (given $\theta$) and we use a conjugate prior
distribution, another (easier) method can be used to determine the
predictive distribution.

Using Bayes Theorem, the posterior density for $\theta$ given
$\xvec$ and $y$ is
\begin{align*}
\pi(\theta|\xvec,y)&=\frac{\pi(\theta)f(\xvec,y|\theta)}{f(\xvec,y)}\\
&=\frac{\pi(\theta)f(\xvec|\theta)f(y|\theta)}{f(\xvec)f(y|\xvec)}
\quad\text{\hfill since $\Xvec$ and $Y$ are independent given $\theta$}\\
&=\frac{\pi(\theta|\xvec)\,f(y|\theta)}{f(y|\xvec)}.
\intertext{Rearranging, we obtain}
f(y|\xvec)&=\frac{f(y|\theta)\pi(\theta|\xvec)}{\pi(\theta|\xvec,y)}.
\end{align*}
This is known as Candidate's formula.  The right-hand-side of this
equation looks as if it depends on $\theta$ but, in fact, any terms
in $\theta$ will be cancelled between the numerator and
denominator.

\begin{example}~\\
  Rework Example~\ref{ex:botulism2} using Candidate's formula to
  determine the number of cases in 2006.

\begin{soln}
  Let $Y$ denote the number of cases in 2006.  We know that
  $\theta|\xvec\sim Ga(8,9)$ and $Y|\theta\sim Po(\theta)$. Using
  Example~\ref{ex:poisgam} we obtain
\[
\theta|\xvec,y\sim Ga(8+y,10).
\]
Therefore the predictive probability function of $Y$ is, for
$y=0,1,\ldots$
\begin{align*}
f(y|\xvec)&=\frac{f(y|\theta)\,\pi(\theta|\xvec)}{\pi(\theta|\xvec,y)}\\
\noalign{\medskip}
&=\bigfrac{\dfrac{\theta^y\,e^{-\theta}}{y!}\times 
\dfrac{9^8\theta^7e^{-9\theta}}{\Gamma(8)}}
{\dfrac{10^{8+y}\theta^{7+y}e^{-10\theta}}{\Gamma(8+y)}}\\ 
\end{align*}
\begin{align*}
\phantom{f(y|\xvec)}
&=\dfrac{\Gamma(8+y)}{y!\Gamma(8)}\times\dfrac{9^8}{10^{8+y}}\\
&=\dfrac{(y+7)!}{y!7!}\times 0.9^8\times 0.1^y\\
&=\binom{y+7}{7}\times 0.9^8\times 0.1^y.
\end{align*}
\end{soln}
\end{example}
%%Hack for vspace
\newpage

\section{Mixture prior distributions}
\label{sec:mixture}
Sometimes prior beliefs cannot be adequately represented by a
simple distribution, for example, a normal distribution or a beta
distribution. In such cases, mixtures of distributions can be useful.

\begin{example}~\\
  Investigations into infants suffering from severe {\it idiopathic
    respiratory distress syndrome} have shown that whether the infant
  survives may be related to their weight at birth. Suppose that you
  are interested in developing a prior distribution for the mean birth
  weight~$\mu$ of such infants. You might have a normal
  $N(2.3,0.52^2)$ prior distribution for the mean birth weight (in kg)
  of infants who survive and a normal $N(1.7,0.66^2)$ prior
  distribution for infants who die. If you believe that the proportion
  of infants that survive is 0.6, what is your prior distribution of
  birth weights of infants suffering from this syndrome?

\begin{soln}
  Let $T=1,2$ denote whether the infant survives or dies. Then the
  information above tells us
\[
\mu|T=1\sim N(2.3,0.52^2)\quad\text{and}\quad \mu|T=2\sim N(1.7,0.66^2).
\]
In terms of density functions, we have
\[
\pi(\mu|T=1)=\phi\left(\mu|2.3,0.52^2\right)
\quad\text{and}\quad 
\pi(\mu|T=2)=\phi\left(\mu|1.7,0.66^2\right),
\]
where $\phi(\cdot|a,b^2)$ is the normal $N(a,b^2)$ density function.

The prior distribution of birth weights of infants suffering from this
syndrome is the (marginal) distribution of~$\mu$. Using the Law of
Total Probability, the marginal density of $\mu$ is
\begin{align*}
  \pi(\mu)&=Pr(T=1)\times\pi(\mu|T=1)+Pr(T=2)\times\pi(\mu|T=2)\\
  &=0.6\,\phi\left(\mu|2.3,0.52^2\right)+0.4\,\phi\left(\mu|1.7,0.66^2\right).
\end{align*}
We write this as
\[
\mu\sim 0.6\,N(2.3,0.52^2)+0.4\,N(1.7,0.66^2).
\]
\end{soln}
%%Hack for vspace
\vspace{-4cm}
\newpage
This prior distribution is a mixture of two normal distributions.
Figure~\ref{fig:mixtureplot} shows the overall (mixture) prior
distribution $\pi(\mu)$ and the ``component'' distributions describing
prior beliefs about the mean weights of those who survive and those
who die. Notice that, in this example, although the mixture
distribution is a combination of two distributions, each with one
mode, this mixture distribution has only one mode. Also, although the
component distributions are symmetric, the mixture distribution is not
symmetric.

\begin{figure}[h]
\begin{center}
\includegraphics[angle=270,width=12cm]{mixtureplot}
\caption{Plot of the mixture density (solid) with its component
densities (survive -- dashed; die -- dotted)}
\label{fig:mixtureplot}
\end{center}
\end{figure}
\end{example}

\begin{definition}~\\
A {\it mixture} of the distributions $\pi_i(\theta)$ with {\it
weights} $p_i$ ($i=1,2,\ldots,m$) has probability (density) function
\begin{equation}
\label{eq:mixture}
\pi(\theta)=\sum_{i=1}^m p_i\pi_i(\theta).
\end{equation}

Figure~\ref{fig:mixtureplot2} contains a plot of two quite different
mixture distributions. One mixture distribution has a single mode and
the other has two modes. In general, a mixture distribution whose $m$
component distributions each have a single mode will have at most $m$
modes.

\begin{figure}[h]
\begin{center}
\includegraphics[angle=270,width=12cm]{mixtureplot2}
\caption{Plot of two mixture densities: solid is
$0.6N(1,1)+0.4N(2,1)$; dashed is $0.9 Exp(1)+0.1N(2,0.25^2)$}
\label{fig:mixtureplot2}
\end{center}
\end{figure}

In order for a mixture distribution to be proper, we must have
\begin{align*}
1&=\int_\Theta \pi(\theta)\,d\theta \\
&=\int_\Theta \sum_{i=1}^m p_i\pi_i(\theta)\,d\theta \\
&=\sum_{i=1}^m p_i\int_\Theta \pi_i(\theta)\,d\theta \\
&=\sum_{i=1}^m p_i, 
\end{align*}
that is, the sum of the weights must be one.
\end{definition}

We can calculate the mean and variance of a mixture distribution as
follows. We will assume, for simplicity, that $\theta$ is a
scalar. Let $E_i(\theta)$ and $Var_i(\theta)$ be the mean and variance
of the distribution for $\theta$ in component~$i$, that is,
\[
E_i(\theta)=\int_\Theta \theta\,\pi_i(\theta)\,d\theta
\quad\text{and}\quad
Var_i(\theta)=\int_\Theta \left\{\theta-E_i(\theta)\right\}^2
\pi_i(\theta)\,d\theta.
\]
It can be shown that the mean of the mixture distribution is
\begin{align}
\label{eq:mixmean}
E(\theta)&=\sum_{i=1}^m p_iE_i(\theta).\\
\intertext{We also have}
E(\theta^2)&=\sum_{i=1}^m p_iE_i(\theta^2)\notag\\
\label{eq:mixexp2}
&=\sum_{i=1}^m p_i\left\{Var_i(\theta)+E_i(\theta)^2\right\}\\
\intertext{from which we can calculate the variance of the mixture
distribution using}
Var(\theta)&=E(\theta^2)-E(\theta)^2\notag.
\end{align}

Combining a mixture prior distribution with data $\xvec$ using Bayes
Theorem produces the posterior density
\begin{align}
\pi(\theta|\xvec)&=\frac{\pi(\theta)\,f(\xvec|\theta)}{f(\xvec)}\notag\\
\label{eq:mixpost}
&=\sum_{i=1}^m\frac{p_i\pi_i(\theta)\,f(\xvec|\theta)}{f(\xvec)}\\
\intertext{where $f(\xvec)$ is a constant with respect to~$\theta$.
  Now if the prior density were $\pi_i(\theta)$ (instead of the
  mixture distribution), using Bayes Theorem, the posterior density
  would be}
\pi_i(\theta|\xvec)&=\frac{\pi_i(\theta)\,f(\xvec|\theta)}{f_i(\xvec)}\notag\\
\intertext{where $f_i(\xvec)$, $i=1,2,\ldots,m$ are constants with
  respect to~$\theta$, that is $\pi_i(\theta)\,f(\xvec|\theta)=f_i(\xvec)\,\pi_i(\theta|\xvec)$.
 Substituting this in to $\eqref{eq:mixpost}$
  gives}
\pi(\theta|\xvec)&=\sum_{i=1}^m\frac{p_if_i(\xvec)}{f(\xvec)}\,\pi_i(\theta|\xvec).\notag\\
\intertext{Thus the posterior distribution is a mixture distribution
  of component distributions $\pi_i(\theta|\xvec)$ with weights
  $p_i^*=p_if_i(\xvec)/f(\xvec)$. Now}
\sum_{i=1}^mp_i^*&=1\quad\Rightarrow\quad
\sum_{i=1}^m\frac{p_if_i(\xvec)}{f(\xvec)}=1\quad\Rightarrow\quad
f(\xvec)=\sum_{i=1}^m p_if_i(\xvec)\notag\\
\intertext{and so} p_i^*&=\bigfrac{p_if_i(\xvec)}{\sum_{j=1}^m
  p_jf_j(\xvec)}, \qquad i=1,2,\ldots,m.\notag
\end{align}

Hence, combining data $\xvec$ with a mixture prior distribution
$(p_i,\pi_i(\theta))$ produces a posterior mixture distribution
$(p_i^*,\pi_i(\theta|\xvec))$. The effect of introducing the data
is to ``update'' the mixture weights ($p_i\to p_i^*$) and the
component distributions ($\pi_i(\theta)\to\pi_i(\theta|\xvec))$.

\begin{example}~\\
Suppose we have a random sample of size 20 from an exponential
\label{ex:gammix}
distribution, that is, $X_i|\theta\sim Exp(\theta)$, $i=1,2,\ldots,20$
(independent). Also suppose that the prior distribution for $\theta$
is the mixture distribution
\begin{equation*}
\theta\sim 0.6\,Ga(5,10)+0.4\,Ga(15,10),
\end{equation*}
as shown in Figure~\ref{fig:mixtureplot3}.  Here the component
distributions are $\pi_1(\theta)=Ga(5,10)$ and
$\pi_2(\theta)=Ga(15,10)$, with weights $p_1=0.6$ and $p_2=0.4$.

\begin{figure}[h!]
\begin{center}
\includegraphics[angle=270,width=12cm]{mixtureplot3}
\caption{Plot of the mixture prior density}
\label{fig:mixtureplot3}
\end{center}
\end{figure}

Using \eqref{eq:mixmean}, the prior mean is
\begin{align*}
E(\theta)&=
0.6\times\frac{5}{10}+0.4\times\frac{15}{10}=0.9\\
\intertext{and, using \eqref{eq:mixexp2}, the prior second moment for
$\theta$ is}
E(\theta^2)&=
0.6\times\frac{5\times 6}{10^2}+0.4\times\frac{15\times
16}{10^2}=1.14\\
\intertext{from which we calculate the prior variance as}
Var(\theta)&=E(\theta^2)-E(\theta)^2=1.14-0.9^2=0.33\\
\intertext{and prior standard deviation as}
SD(\theta)&=\sqrt{Var(\theta)}=\sqrt{0.33}=0.574.
\end{align*}

We have already seen that combining a random sample of size~20 from an
exponential distribution with a $Ga(g,h)$ prior distribution results
in a $Ga(g+20,h+20\bar x)$ posterior distribution. Therefore, the
(overall) posterior distribution will be a mixture distribution with
component distributions
\[
\pi_1(\theta|\xvec)=Ga(25,10+20\bar x)\quad\text{and}\quad
\pi_2(\theta|\xvec)=Ga(35,10+20\bar x). 
\]
We now calculate new values for the weights $p_1^*$ and
$p_2^*=1-p_1^*$, which will depend on both prior information and the
data. We have
\begin{align*}
p_1^*
&=\frac{0.6f_1(\xvec)}{0.6f_1(\xvec)+0.4f_2(\xvec)}\\
\intertext{from which}
\left(p_1^*\right)^{-1}-1
&=\frac{0.4f_2(\xvec)}{0.6f_1(\xvec)}.
\end{align*}
In general, the functions 
\[
f_i(\xvec)=\int_\Theta\pi_i(\theta)\,f(\xvec|\theta)\,d\theta
\]
are potentially complicated integrals (solved either analytically or
numerically). However, as with Candidates formula, these calculations
become much simpler when we have a conjugate prior distribution:
rewriting Bayes Theorem, we obtain
\[
f(\xvec)=\frac{\pi(\theta)\,f(\xvec|\theta)}{\pi(\theta|\xvec)}
\]
and so when the prior and posterior densities have a simple form (as
they do when using a conjugate prior), it is straightforward to
determine $f(\xvec)$ using algebra rather than having to use calculus.

In this example we know that the gamma distribution is the conjugate
prior distribution: using a random sample of size $n$ with mean $\bar
x$ and a $Ga(g,h)$ prior distribution gives a $Ga(g+n,h+n\bar x)$
posterior distribution, and so
\begin{align*}
f(\xvec)
&=\frac{\pi(\theta)\,f(\xvec|\theta)}{\pi(\theta|\xvec)}\\ 
\noalign{\medskip}
&=\bigfrac{\frac{h^g\theta^{g-1}e^{-h\theta}}{\Gamma(g)}\times
\theta^ne^{-n\bar x\theta}}
{\frac{(h+n\bar x)^{g+n}\theta^{g+n-1}
e^{-(h+n\bar x)\theta}}{\Gamma(g+n)}}\\ 
\noalign{\medskip}
&=\frac{h^g\,\Gamma(g+n)}{\Gamma(g)(h+n\bar x)^{g+n}}.
\end{align*}
Therefore
\begin{align*}
\left(p_1^*\right)^{-1}-1
&=\left.\frac{0.4\times 10^{15}\,\Gamma(35)}{\Gamma(15)(10+20\bar x)^{35}}
\right/\frac{0.6\times 10^5\,\Gamma(25)}{\Gamma(5)(10+20\bar x)^{25}} \\
&=\frac{2\Gamma(35)\Gamma(5)}{3\Gamma(25)\Gamma(15)(1+2\bar x)^{10}} \\
&=\frac{611320}{7(1+2\bar x)^{10}}
\end{align*}
and so
\[
p_1^*=\frac{1}{1+\dfrac{611320}{7(1+2\bar x)^{10}}},
\quad\quad p_2^*=1-p_1^*.
\]

Hence the posterior distribution is the mixture distribution
\[
\frac{1}{1+\dfrac{611320}{7(1+2\bar x)^{10}}}
\times Ga(25,10+20\bar x) +
\left(1-\frac{1}{1+\dfrac{611320}{7(1+2\bar x)^{10}}}\right)
\times Ga(35,10+20\bar x).
\]
Recall that the most likely value of $\theta$ from the data alone, the
likelihood mode, is $1/\bar x$. Therefore, large values of $\bar x$
indicate that $\theta$ is small and {\it vice versa}. With this in
mind, it is not surprising that the weight $p_1^*$ (of the component
distribution with the smallest mean) is increasing in $\bar x$, and
$p_1^*\to 1$ as $\bar x\to\infty$. Using \eqref{eq:mixmean}, the
posterior mean is
\begin{align*}
E(\theta|\xvec) 
&=\frac{1}{1+\dfrac{611320}{7(1+2\bar x)^{10}}}
\times\dfrac{25}{10+20\bar x} +
\left(1-\frac{1}{1+\dfrac{611320}{7(1+2\bar x)^{10}}}\right)
\times\frac{35}{10+20\bar x} \\
&=\cdots\\
&=\frac{1}{2(1+2\bar x)}
\left\{7-\frac{2}{1+\dfrac{611320}{7(1+2\bar x)^{10}}}\right\}.
\end{align*}
The posterior standard deviation can be calculated using
\eqref{eq:mixmean} and \eqref{eq:mixexp2}.

Table~\ref{tab:mixture} shows the posterior distributions which result
when various sample means $\bar x$ are observed together with the
posterior mean and the posterior standard deviation. Graphs of these
posterior distributions, together with the prior distribution, are
given in Figure~\ref{fig:mixtureplot4}. When considering the effect on
beliefs of observing the sample mean $\bar x$, it is important to
remember that large values of $\bar x$ indicate that $\theta$ is small
and {\it vice versa}. Plots of the posterior mean against the sample
mean reveal that the posterior mean lies between the prior mean and
the likelihood mode only for $\bar x\in(0,0.70)\cup(1.12,\infty)$. Note that observing the data has
focussed our beliefs about $\theta$ in the sense that the posterior
standard deviation is less than the prior standard deviation -- and
considerably less in some cases.

\newpage

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
$\bar x$ & $\hat\theta=1/\bar x$ & Posterior mixture distribution &
$E(\theta|\xvec)$ & $SD(\theta|\xvec)$\\
\hline
4 & 0.25 & $0.99997\,Ga(25,90)+0.00003\,Ga(35,90)$ & 0.278 & 0.056\\
2 & 0.5 & $0.9911\,Ga(25,50)+0.0089\,Ga(35,50)$ & 0.502 & 0.102\\
1.2 & 0.8 & $0.7027\,Ga(25,34)+0.2973\,Ga(35,34)$ & 0.823 & 0.206\\
1 & 1.0 & $0.4034\,Ga(25,30)+0.5966\,Ga(35,30)$ & 1.032 & 0.247\\
0.8 & 1.25 & $0.1392\,Ga(25,26)+0.8608\,Ga(35,26)$ & 1.293 & 0.260\\
0.5 & 2.0 & $0.0116\,Ga(25,20)+0.9884\,Ga(35,20)$ & 1.744 & 0.300\\
\hline
\end{tabular}
\caption{Posterior distributions (with summaries) for various sample
means $\bar x$} 
\label{tab:mixture}
\end{center}
\end{table}

\begin{figure}[h!]
\begin{center}
\includegraphics[angle=270,width=12cm]{mixtureplot4}
\caption{Plot of the prior distribution and various posterior
distributions} 
\label{fig:mixtureplot4}
\end{center}
\end{figure}
\end{example}

%%Hack for vspace

\newpage

\section{Learning objectives}
By the end of this chapter, you should be able to:
\begin{itemize}
\item determine the likelihood function using a random
  sample from \textbf{any} distribution
\item combine this likelihood function with \textbf{any}
  prior distribution to obtain the posterior distribution
\item name the posterior distribution if it is a ``standard''
  distribution listed in these notes or on the exam paper -- this list
  may well include distributions that are standard within the subject
  but which you have not met before. If the posterior distribution is
  not a ``standard'' distribution then it is okay just to give its
  density (or probability function) up to a constant.
\item do all the above for a particular data set or for a general case
  with random sample $x_1,\ldots,x_n$
\item describe the different levels of prior information;
  determine and use conjugate priors and vague priors
\item determine the asymptotic posterior distribution
\item determine the predictive distribution, particularly when having
  a random sample from any distribution and a conjugate prior via
  Candidate's formula
\item describe and calculate the confidence intervals, HDIs and
  prediction intervals
\item determine posterior distributions when the prior is a mixture of
  conjugate distributions
\end{itemize}