\documentclass[12pt,a4paper]{article} \usepackage{amsthm,amsmath,amsfonts} \usepackage{latexsym} \usepackage{fancyhdr} \usepackage{enumerate} %\usepackage{marginnote} \usepackage{graphicx} \usepackage{skull} \usepackage{hyperref} \usepackage{fancybox} \usepackage{xspace} %%From web for code highlighting:: \usepackage{listings} \usepackage{cmbright} \addtolength{\parskip}{1ex} \baselineskip = 6mm %\parskip = 2mm \textheight 23.6cm \textwidth 16.6cm \oddsidemargin 0.0cm \evensidemargin 0.0cm \topmargin -1.4cm %\headsep -.60in \setlength{\parindent}{0cm} \pagestyle{fancy} \lhead{} \chead{} \rhead{} \lfoot{} \cfoot{\thepage} \rfoot{} \rhead{} \newcommand{\pr}[2][]{\ensuremath{\mathrm{Pr}_{#1}\!\left(#2\right)}} \newenvironment{myhighlight}{\textit}{} \renewenvironment{emph}{\textbf}{} \newcommand{\Exp}[2][]{\ensuremath{\mathrm{Exp}\left(#2\right)}} \newcommand{\Bin}[2][]{\ensuremath{\mathrm{Bin}\left(#2\right)}} \newcommand{\Po}[2][]{\ensuremath{\mathrm{Po}\left(#2\right)}} \newcommand{\cov}[2][]{\ensuremath{\mathrm{Cov}\left(#2\right)}} \newcommand{\cor}[2][]{\ensuremath{\mathrm{Cor}\left(#2\right)}} \newcommand{\diag}[2][]{\ensuremath{\mathrm{diag}\left(#2\right)}} \newcommand{\given}{\ensuremath{\ |\ }} \newcommand{\stationary}{\ensuremath{\vec{\pi}}} \newcommand{\identity}{\ensuremath{\vec{I}}} \renewcommand{\vec}[1]{\underline{#1}} \newcommand{\remove}[1]{} %\newcommand{\cover}[2]{\hi{#1}} \newcommand{\cover}[2]{#2} \newcommand{\hi}[2][]{\textbf{#2}} \newcommand{\bigfont}{\fontsize{1.45cm}{1.5cm}\selectfont} \newcommand{\announce}{\setbeamercolor{structure}{fg=green!60!black}} \newcommand{\sol}{\setbeamercolor{structure}{fg=purple!80!black}} \newcommand{\normal}{\setbeamercolor{structure}{fg=blue!60!green}} \newcommand{\vsp}{\vspace{0.3cm}} \newcommand{\vspl}{\vspace{0.15cm}} \newcommand{\xvec}{\boldsymbol{x}} \newcommand{\yvec}{\vec{y}} \newcommand{\Xvec}{\boldsymbol{X}} \newcommand{\Yvec}{\vec{Y}} \newcommand{\Evec}{\vec{E}} \newcommand{\tvec}{\bodsymbol{t}} \newcommand{\Tvec}{\bodsymbol{T}} \newcommand{\mutau}{\begin{pmatrix}\mu\\\tau\end{pmatrix}} \newcommand{\thetavec}{\vec{\theta}} \newcommand{\thetavechat}{\vec{\hat\theta}} \newcommand{\thetahat}{\hat\theta} \newcommand{\xvece}{\boldsymbol{x}=(x_1,x_2,\ldots,x_n)^T} \newcommand{\Xvece}{\boldsymbol{X}=(X_1,X_2,\ldots,X_n)^T} \newcommand{\thetavece}{\vec{\theta}=(\theta_1,\theta_2,\ldots,\theta_p)^T} \newcommand{\exxvec}{$x_1,x_2,\ldots,x_n$\xspace} \newcommand{\exxvecdot}{$x_1,x_2,\ldots,x_n$.\ } \newcommand{\exXvec}{$X_1,X_2,\ldots,X_n$\ } \newcommand{\muvec}{\vec{\mu}} \newtheoremstyle{exmp}% hnamei {3pt}% hSpace abovei {3pt}% hSpace belowi {}% hBody fonti {}% hIndent amounti1 {\bf}% hTheorem head fonti {:}% hPunctuation after theorem headi {.5em}% hSpace after theorem headi2 {}% hTheorem head spec (can be left empty, meaning `normal')i \theoremstyle{exmp} \newtheorem{example}{\textcolor{purple}{\textbf{Example}}}[section] %\def \var{{\textrm{Var}}} \newcommand{\E}[2][]{\ensuremath{\mathrm{E}\left[#2\right]}} \def \R{\mathbb{R}} \def \Z{\mathbb{Z}} \def \proof{\noindent{\bf Proof: }} \def \eproof{\qed\vskip .3cm} \def \union{\cup} \def \intersect{\cap} \newcommand{\ex}[2][]{\ensuremath{\mathrm{E}_{#1}\!\left[#2\right]}} \newcommand{\var}[2][]{\ensuremath{\mathrm{Var}_{#1}\left(#2\right)}} %\renewcommand{\det}[1][]{\ensuremath{\mathrm{det}\left(#1\right)}} \newtheorem{defn}{Definition} %\newtheorem*{example}{Example} \newtheorem{result}{Result}[section] \newtheorem{algorithm}{Algorithm}[section] \usepackage{Sweave} \usepackage[usenames, dvipsnames]{xcolor} \begin{document} %%From web for code/highlighting \definecolor{mygreen}{rgb}{0,0.6,0} \definecolor{mygray}{rgb}{0.5,0.5,0.5} \definecolor{mymauve}{rgb}{0.58,0,0.82} \definecolor{light-gray}{gray}{0.95} \definecolor{LEEgreen}{rgb}{0.04,0.667,0.137} \lstset{ backgroundcolor=\color{light-gray}, % choose the background color; you must add \usepackage{color} or \usepackage{xcolor}; should come as last argument basicstyle=\footnotesize, % the size of the fonts that are used for the code breakatwhitespace=false, % sets if automatic breaks should only happen at whitespace breaklines=true, % sets automatic line breaking captionpos=b, % sets the caption-position to bottom commentstyle=\color{mygreen}, % comment style deletekeywords={...}, % if you want to delete keywords from the given language escapeinside={\%*}{*)}, % if you want to add LaTeX within your code extendedchars=false, % lets you use non-ASCII characters; for 8-bits encodings only, does not work with UTF-8 frame=single, % adds a frame around the code keepspaces=true, % keeps spaces in text, useful for keeping indentation of code (possibly needs columns=flexible) keywordstyle=\color{blue}, % keyword style language=R, % the language of the code morekeywords={*,...}, % if you want to add more keywords to the set numbers=left, % where to put the line-numbers; possible values are (none, left, right) numbersep=5pt, % how far the line-numbers are from the code numberstyle=\tiny\color{mygray}, % the style that is used for the line-numbers rulecolor=\color{black}, % if not set, the frame-color may be changed on line-breaks within not-black text (e.g. comments (green here)) showspaces=false, % show spaces everywhere adding particular underscores; it overrides 'showstringspaces' showstringspaces=false, % underline spaces within strings only showtabs=false, % show tabs within strings adding particular underscores stepnumber=1, % the step between two line-numbers. If it's 1, each line will be numbered stringstyle=\color{mymauve}, % string literal style tabsize=2, % sets default tabsize to 2 spaces title=\lstname % show the filename of files included with \lstinputlisting; also try caption instead of title } %\SweaveOpts{concordance=TRUE} \pagenumbering{gobble} \setcounter{page}{-3} % to get contents page labelled (i) %%%Practicals %\section*{Aims} %This introductory practical is split into two parts: %\begin{itemize} %\item In the first hour or so, we will focus on the basics of R, revising some of the key concepts from last year's MAS1802 material %\item In the second hour we will focus on material from this week's lecture -- using R to simulate from standard probability distributions, as well as more complex models including mixture distributions %\end{itemize} %\vspace{0.25cm} %\noindent Please read through all of this handout \textit{very carefully}, and follow the instructions \textit{meticulously}. Whenever %you are typing commands into R, make sure you do so \textit{extremely carefully}. If you get an error message, %try typing the code in again, just to make sure you've not made a spelling mistake or other simple typo. A %common error is to forget that R is case sensitive, so use upper/lower cases exactly when you're told to %do so. %\vspace{0.25cm} %\noindent For all R practicals, make sure you've got your lecture notes handy. %\section*{Getting started} %\noindent Log-in to a PC in the usual way. Now open Rstudio. If you can't remember how to do this, click on the ``Start'' button in the bottom %left-hand-corner and search for ``rstudio''. Once you've opened Rstudio you should see four panels; if there is no text editor (usually the top-left panel), click on the green button with a white plus symbol and open a new R script. Remember -- to execute a line of code in the text editor, hit \texttt{Ctrl}-\texttt{Enter}. %\vspace{0.25cm} %\noindent It is good practice to always use the text editor, rather than type into the console window (bottom- %left panel) directly. This is especially true in future practicals when you may be building large functions; it %is much easier to de-bug and edit code in the text editor than in the console window. %\vspace{0.25cm} %\noindent \textbf{Hint: To commit several lines of code in one go, highlight your block of code in the text editor and hit \texttt{Ctrl}-\texttt{Enter} just once!} %\vspace{0.25cm} %\noindent Now you've got Rstudio open, you're ready to go! Try out \textit{all} of the following exercises. %\newpage \lhead{MAS3902} \rhead{Bayesian Inference} \section*{Bayesian Inference: Assignment 3} \noindent Answers to the starred (\textbf{*}) questions should be submitted to Canvas by no later than \textcolor{red}{\bf 4pm, Friday 22nd March 2024}. Solutions to non-assessed questions will be made available in class and/or in Canvas. \subsection*{Question 1* (30 marks total)} A hepatologist is interested in the levels of the liver enzyme \textit{ornithine carbonyltransferase} in patients suffering from acute viral hepatitis. She collects measurements from a random sample of patients and the logarithm of their enzyme measurements are given in the following table. They are also available in the \texttt{R} datafile \texttt{hepatitis} in the \texttt{nclbayes} package; alternatively, you can download these data from a webpage by typing (but \textbf{not} copy-and-pasting!): {\small \begin{lstlisting} hepatitis = read.table("https://www.mas.ncl.ac.uk/~nlf8/teaching/mas3902/hepatitis.txt") \end{lstlisting} } \vspace{-1cm} \begin{center} \begin{tabular}{|cccccccccc|} \hline 2.64 & 2.51 & 2.20 & 2.53 & 2.02 & 2.47 & 2.75 & 2.77 & 2.91 & 2.45 \\ 2.25 & 1.96 & 2.22 & 2.23 & 1.98 & 2.70 & 2.61 & 2.76 & 2.03 & 2.38 \\ 2.62 & 2.28 & 2.47 & 3.04 & 1.91 & 2.71 & 2.89 & 2.70 & 2.29 & 2.50 \\ \hline \end{tabular} \end{center} Assume that the enzyme measurement varies according to a $N(\mu,1/\tau)$ distribution. After extensiver discussions with an expert, prior beliefs about $\mu$ and $\tau$ are summarised as \[ \mutau\sim NGa(2.6,1,5,0.4). \] \begin{enumerate}[(a)] %\item Use a normal probability (q-q) plot to confirm the suitability % of the normal distribution as a model for the variation in enzyme % measurements. The relevant \texttt{R} commands are \texttt{qqnorm} % and \texttt{qqline}. \item Calculate the expert's prior mean and standard deviation for $\mu$, $\tau$ and $\sigma=1/\sqrt\tau$. \item Determine the posterior distribution for $(\mu,\tau)^T$ after combining the hepatologist's prior beliefs with the data. Calculate the posterior mean and standard deviation of $\mu$, $\tau$ and $\sigma$. \item From the lecture notes, we know that $\mu|\boldsymbol{x}\sim t_{2G}(B, H/GC)$ and $\tau|\boldsymbol{x}\sim Ga(G, H)$ (where expressions for $B$, $C$, $G$ and $H$ are given in the notes). Show that $\sigma|\boldsymbol{x}\sim \textit{Inv-Chi}(G, H)$. \item Plot the (marginal) prior and posterior densities for $\mu$ on the same graph. Construct similar plots for $\tau$ and $\sigma$. Also produce contour plots of the (joint) prior and posterior densities for $(\mu,\tau)^T$ on the same graph. \item Plot 80\%, 90\% and 95\% prior and posterior confidence regions for $(\mu,\tau)^T$ on the same graph. \item Use these plots and your calculations to comment on the main changes in the hepatologist's beliefs about $\mu$, $\tau$ and $\sigma$ after incorporating the data. Include a comment on the prior-to-posterior change in the dependence structure (contour shape) of $(\mu,\tau)$ and on their confidence regions for $(\mu,\tau)$. \item The hepatologist is particularly interested in whether the population mean level $\mu$ is larger than 2.7. Determine the prior and posterior probabilities for $\mu>2.7$. Have the data been informative? \item The hepatologist starts to think about the enzyme level $Y$ in the next patient admitted to the study. Write down the predictive distribution for this future observation, and find the predictive probability that $Y>2.7$. \end{enumerate} \subsection*{Question 2} Suppose you are interested in the aerobic capacity (ml/kg/min) of athletes, measurements on which are assumed Normally distributed with mean $\mu$ and precision $\tau$. \textit{Prior specification:} \begin{itemize} \item [(a)] In a previous study, the aerobic capacity of athletes was assumed Normal with known standard deviation $\sigma=\sqrt{3}$ ml/kg/min; the mean $\mu$ was given a Normal prior such that $\mu \sim N(71,1)$. Assuming the same values for prior mean and variance in the current study, find $b$ and $c$ such that \[ \mu|\tau \sim N\left(b, \frac{1}{c\tau}\right). \] \item [(b)] Unlike the previous study, in the current study we are assuming \textit{both} $\mu$ \textit{and} $\tau$ (and hence $\sigma=1/\sqrt{\tau}$) are unknown. We assume a Gamma prior for the precision $\tau$, that is, \[ \tau \sim Ga(g, h). \] Discussions with a sports scientist suggest that $E(\tau) = Var(\tau) = 2$. Find $g$ and $h$. \item [(c)] Using results from Chapter 2, write down the joint prior for $(\mu, \tau)^{T}$. \end{itemize} \textit{Posterior distribution:} The aerobic capacity of 50 randomly chosen athletes was measured, giving a mean of $\bar{x} = 71.5$ ml/kg/min and standard deviation $s = 2.3$ ml/kg/min. \begin{itemize} \item [(d)] Using results from Chapter 2, find the posterior distribution for $(\mu, \tau)^{T}$. \item [(e)] Complete the following table: {\scriptsize \begin{center} \begin{tabular}{|c|cc|cc|} \hline $\mu$ &\multicolumn{2}{l|}{Prior distribution: $\mu \sim$\textcolor{white}{answer here}} &\multicolumn{2}{l|}{Posterior distribution: $\mu|\boldsymbol{x} \sim$\textcolor{white}{answer here}} \\ & & & & \\ &Mean: \textcolor{white}{Answer} &Variance: \textcolor{white}{Answer} &Mean: \textcolor{white}{Answer} &Variance: \textcolor{white}{Answer}\\ & & & & \\ \hline $\tau$ &\multicolumn{2}{l|}{Prior distribution: $\tau \sim$\textcolor{white}{answer}} &\multicolumn{2}{l|}{Posterior distribution: $\tau|\boldsymbol{x} \sim$ \textcolor{white}{answer goes here}} \\ & & & & \\ &Mean: \textcolor{white}{Answer} &Variance: \textcolor{white}{Answer} &Mean: \textcolor{white}{Answer} &Variance: \textcolor{white}{Answer} \\ & & & & \\ \hline $\sigma$ &\multicolumn{2}{l|}{Prior distribution: $\sigma \sim$\textcolor{white}{answer}} &\multicolumn{2}{l|}{Posterior distribution: $\sigma|\boldsymbol{x}\sim$\textcolor{white}{answer}} \\ & & & & \\ &Mean: \textcolor{white}{Answer} &Variance: \textcolor{white}{Answer} &Mean: \textcolor{white}{Answer} &Variance: \textcolor{white}{Answer} \\ & & & & \\ \hline \end{tabular} \end{center} } \item [(f)] Sketch the prior and posterior for $\mu$ on the same graph. \item [(g)] Which of the plots in Figure 1 (a) -- left or right -- correspond to the prior (dashed) and posterior (bold) distributions for $(\mu,\tau)^{T}$? Explain your answer. \end{itemize} \textit{Confidence intervals and regions:} \begin{itemize} \item[(h)] Find the 95\% highest density intervals for $\mu|\boldsymbol{x}$. You might find the following output from R useful: \begin{verbatim} > qt(0.975,54) [1] 2.004879 \end{verbatim} \item [(i)] Why are highest density intervals for $\tau$ and $\sigma$ not as straightforward to obtain? Use a Normal approximation to obtain approximate 95\% confidence intervals for $\tau|\boldsymbol{x}$ and $\sigma|\boldsymbol{x}$. \item [(j)] Describe what the plot in Figure 1(b) shows. \end{itemize} \textit{Prediction:} Suppose we are interested in the aerobic capacity of the next athlete to join the study. \begin{itemize} \item [(k)] Using results from Chapter 2, find the predictive distribution for the aerobic capacity of this new athlete. \item [(l)] Find the 95\% highest density interval associated with the predictive distribution in (i). \end{itemize} \begin{figure}[!h] \centering \includegraphics[width=8cm,height=6.5cm]{extraQ.pdf} \includegraphics[width=8cm, height=6.5cm]{cregion.pdf} \caption{(a) Prior (dashed) and posterior (solid) densities (left); (b) contours from a bivariate density (right).} \end{figure} \newpage \section*{Solutions} \subsection*{Question 2} \begin{itemize} \item [(a)] From the previous study, we have $X_{i}|\mu \sim N(\mu, 3)$ (so $\tau=1/3$); we also have that $\mu \sim N(71, 1)$. In \textit{our} study, we assume that $\mu|\tau \sim N(b, 1/c\tau)$. \begin{itemize} \item Keeping the prior means the same gives $b=71$ \item Keeping the prior variances the same gives $1/c\tau = 1$, that is, \[ \frac{1}{c\left(\frac{1}{3}\right)} = 1, \] and so $c=3$. \end{itemize} \item [(b)] We assume $\tau \sim Ga(g, h)$. \begin{itemize} \item Using $E(\tau) = g/h= 2$ gives $g=2h$. \item Using $Var(\tau) = g/h^{2} = 2$ gives \[ \frac{2h}{h^{2}} = 2 \longrightarrow h=1; \] therefore $g=2$. \end{itemize} \item [(c)] So, from (a) and (b) we have that \[ \left(\begin{array}{c} \mu \\ \tau \end{array} \right) \sim NGa(71, 3, 2, 1). \] \item [(d)] From chapter 2, we know that \[ \left(\begin{array}{c} \mu \\ \tau \end{array} \right)|\boldsymbol{x} \sim NGa(B, C, G, H), \] with $B$, $C$, $G$, and $H$ as defined on page 32 of the lecture notes. With $b=71, c=3, g=2, h=1$ and $n=50, \bar{x}=71.5, s=2.3$, we have \[ B = 71.4717, C = 53, G = 27, H = 133.6038. \] Therefore, \[ \left(\begin{array}{c} \mu \\ \tau \end{array} \right)|\boldsymbol{x} \sim NGa(71.4717, 53, 27, 133.6038), \] \item[(e)] For example, we have $\mu|\boldsymbol{x} \sim t_{2G}(B, H/GC) = t_{54}(71.4717, 0.0934)$, with \[ E(\mu|\boldsymbol{x}) = 71.4717; \quad Var(\mu|\boldsymbol{x}) = \frac{54 \times 0.0934}{52} = 0.09699, \] using the formula sheet. Similarly, for the prior and the other parameters, gives: \begin{center} \begin{tabular}{|c|cc|cc|} \hline $\mu$ &\multicolumn{2}{l|}{Prior distribution: $\mu \sim t_{4}(71, 0.1667)$} &\multicolumn{2}{l|}{Posterior distribution: $\mu|\boldsymbol{x} \sim t_{54}(71.4717,0.0934)$} \\ & & & & \\ &Mean: 71 &Variance: 0.3334 &Mean: 71.4717 &Variance: 0.09699\\ & & & & \\ \hline $\tau$ &\multicolumn{2}{l|}{Prior distribution: $\tau \sim Ga(2,1)$} &\multicolumn{2}{l|}{Posterior distribution: $\tau|\boldsymbol{x} \sim Ga(27, 133.6038)$} \\ & & & & \\ &Mean: 2 &Variance: 2 &Mean: 0.2021 &Variance: 0.00153 \\ & & & & \\ \hline $\sigma$ &\multicolumn{2}{l|}{Prior distribution: $\sigma \sim InvChi(2,1)$} &\multicolumn{2}{l|}{Posterior distribution: $\sigma|\boldsymbol{x}\sim InvChi(27,133.6038)$} \\ & & & & \\ &Mean: 0.8862 &Variance: 0.2146 &Mean: 2.2560&Variance: 0.0491 \\ & & & & \\ \hline \end{tabular} \end{center} \item [(g)] Must be the left-hand plot -- the modes are at the incorrect place for $\mu$ in both prior and posterior in the right-hand plot. \item [(h)] Using the fact that $\mu|\boldsymbol{x} \sim t_{54}(71.4717, 0.0934)$, we can find \[ 71.4717 \pm 2.00489\sqrt{0.0934} \longrightarrow (70.859, 72.8442). \] \item [(i)] Not usually symmetric! For $\tau|\boldsymbol{x} \sim Ga(27, 133.6038)$, we have $E(\tau|\boldsymbol{x}) = 0.2021$ and $Var(\tau|\boldsymbol{x}) = 0.00153$. Therefore, approximately, we have \[ \tau|\boldsymbol{x} \sim N(0.2021, 0.00153), \] giving $0.2021 \pm 1.96 \times \sqrt{0.00153} \longrightarrow (0.1254, 0.2788)$. Similarly for $\sigma|\boldsymbol{x}$ we have $(1.8217, 2.6903)$. \item[(j)] $80\%$, $90\%$ and $95\%$ credible regions for $(\mu, \tau)^{T}$. \item [(k)] $Y|\boldsymbol{x}\sim t_{54}(71.4717, 5.0417)$. \item [(l)] $71.4717 \pm 2.004879 \times \sqrt{5.0417} \longrightarrow (66.9699, 75.9734)\text{ ml/kg/min.}$. \end{itemize} \end{document} \subsection*{Question 1} \begin{enumerate}[(a)] \item If $\begin{pmatrix}\mu\\\tau\end{pmatrix}\sim NGa(b,c,g,h)$ then the univariate marginal distributions are $\mu\sim t_{2g}\{b,h/(gc)\}$, $\tau\sim Ga(g,h)$ and $\sigma\sim\text{Inv-Chi}(g,h)$. Therefore \textcolor{red}{([0.5 marks] for each of mean and SD, so [3 marks] in total)}: \begin{align*} E(\mu)&=b=2.6,\qquad SD(\mu)=\sqrt{\frac{2gh}{gc(2g-2)}}=\sqrt{\frac{h}{c(g-1)}}=0.3162. \qquad \intertext{Also} E(\tau)&=\frac{g}{h}=12.5,\qquad SD(\tau)=\frac{\sqrt{g}}{h}=5.5902 \qquad \intertext{and} E(\sigma)&=\dfrac{\sqrt{h}\,\Gamma\left(g-\frac{1}{2}\right)}{\Gamma(g)}=0.3065\\ SD(\sigma)&=\sqrt{\frac{h}{g-1}-E(\sigma)^2}=0.0777.\qquad \end{align*} \item The posterior distribution is $(\mu,\tau)^T|\xvec\sim NGa(B,C,G,H)$ where \textcolor{red}{([0.5 marks] for each of $B$, $C$, $G$ and $H$, so [2 marks] in total)}: \begin{align*} B&=\frac{bc+n\bar x}{c+n}=2.46387, C=c+n=31,\\ G&=g+\frac{n}{2}=20, H=h+\frac{cn(\bar x-b)^2}{2(c+n)}+\frac{ns^2}{2}=1.76697. \end{align*} Therefore \textcolor{red}{([0.5 marks] for each of mean and SD, so [3 marks] in total)}: \begin{align*} E(\mu|\xvec)&=B=2.46387, SD(\mu|\xvec)=\sqrt{\frac{H}{C(G-1)}}=0.05476. \intertext{Also} E(\tau|\xvec)&=\frac{G}{H}=11.3252, SD(\tau|\xvec)=\frac{\sqrt{G}}{H}=2.5324 \intertext{and} E(\sigma|\xvec)&=\dfrac{\sqrt{H}\,\Gamma\left(G-\frac{1}{2}\right)}{\Gamma(G)} =0.3029,\qquad SD(\sigma|\xvec)=\sqrt{\frac{H}{G-1}-E(\sigma|\xvec)^2}=0.0349. \end{align*} \item For ease of notation, use $X$ and $Y$ for $\tau$ and $\sigma$ respectively. Then something like: \[ F_{Y}(y) = Pr(Y\leq y) = Pr\left(\frac{1}{\sqrt{X}} \leq y\right) = Pr\left(X \geq \frac{1}{y^{2}}\right) = 1-F_{X}\left(\frac{1}{y^{2}}\right). \] Then differentiate to get the PDF: \[ f_{Y}(y) = \frac{d}{dy}F_{Y}(y) = \frac{d}{dy}\left(1-F_{X}\left(\frac{1}{y^2}\right)\right) = -f_{X}\left(\frac{1}{y^2}\right) \times -2 y^{-3} = 2 f_{X}\left(\frac{1}{y^2}\right) \times y^{-3}. \] We know that $X\sim Ga(G, H)$, and so \[ f_{X}\left(\frac{1}{y^2}\right) = \frac{H^{G} (y^{-2})^{G-1} e^{-Hy^{-2}}}{\Gamma(G)}; \] therefore \begin{eqnarray*} f_{Y}(y) &= &2 \times \frac{H^{G} (y^{-2})^{G-1} e^{-Hy^{-2}}}{\Gamma(G)} \times y^{-3} \\ & & \\ &=& \frac{2H^{G}y^{-2G+2-3} e^{-H/y^2}}{\Gamma(G)} \\ & & \\ &=& \frac{2H^{G}y^{-2G-1} e^{-H/y^2}}{\Gamma(G)}, \end{eqnarray*} which, by comparison with the formula sheet, is the PDF for an $\textit{Inv-Chi}(G, H)$ distribution. \textcolor{red}{[6 marks]; students might choose to use the transformation method formula instead -- this is fine.} \item Plots of the (marginal) prior and posterior densities for $\mu$, $\tau$ and $\sigma$ are given in Figure~\ref{fig:proj1-post}. The contour plot of the (joint) prior and posterior densities for $(\mu,\tau)$ is shown in Figure~\ref{fig:proj1-post-joint}. \textcolor{red}{[4 marks]; plots produced easily in R.} \begin{figure}[t!] \includegraphics[width=0.5\textwidth]{proj1-post-mu}\\ \includegraphics[width=0.5\textwidth]{proj1-post-tau} \includegraphics[width=0.5\textwidth]{proj1-post-sigma} \caption{Prior (dashed) and posterior (solid) densities for $\mu$, $\tau$ and $\sigma$} \label{fig:proj1-post} \end{figure} \begin{figure}[t] \centering \includegraphics[width=0.69\textwidth]{proj1-post-joint} \caption{Contour plot of the prior (dashed) and posterior (solid) densities for $(\mu,\tau)^T$} \label{fig:proj1-post-joint} \end{figure} \item A plot of the 80\%, 90\% and 95\% prior and posterior confidence regions for $(\mu,\tau)^T$ is given in Figure~\ref{fig:proj1-conf}. \textcolor{red}{[2 marks]; easily produced in R, using the \texttt{nclbayes} package.} \begin{figure}[t!] \centering \includegraphics[width=0.69\textwidth]{proj1-conf} \caption{95\%, 90\% and 80\% prior (dashed) and posterior (solid) confidence regions for $(\mu,\tau)^T$} \label{fig:proj1-conf} \end{figure} \item Changes in beliefs: \textcolor{red}{[5 marks] maximum}; \item[] $\mu$: its mean, mode and standard deviation have decreased \item[] $\tau$: its mean has decreased (slightly); its mode has increased (slightly); its standard deviation has more than halved \item[] $\sigma$: its mode has increased slightly; its standard deviation has decreased \item[] $(\mu,\tau)$: the contours are less ``triangular'' and more ``elliptical'' \item[] $(\mu,\tau)$: the confidence regions are much smaller \item The probabilities are $Pr(\mu>2.7)=0.36551$ and $Pr(\mu>2.7|\xvec)=0.000036$, calculated using \texttt{1-pgt(2.7,2*g,b,h/(g*c))} and \texttt{1-pgt(2.7,2*G,B,H/(G*C))}. There is a large difference between these probabilities and so the data have been very informative. \textcolor{red}{[3 marks]; one for each probability, one for the comment.} \item The predictive distribution is given by \[ y|\boldsymbol{x} \sim t_{2G}\left(B, \frac{H(C+1)}{GC}\right) \longrightarrow t_{40}\left(2.46387, \frac{1.76697\times 32}{20\times 31}\right) = t_{40}(2.46387, 0.09120). \quad \text{\textcolor{red}{[1 mark]}} \] Then $Pr(Y>2.7|\xvec) =0.2194$, using \texttt{1-pgt(2.7, 40, 2.46387, 0.09120)}. \textcolor{red}{[1 mark]} \end{enumerate} \newpage The hepatologist starts to think about the enzyme levels in the next sample of $m$ patients. Re-visit question 1 from assignment 1. Recall that we were interested in $\theta$, the probability that particulate matter in the air as a result of a dust storm, results in an environmental health warning being issued. \begin{itemize} \item [(a)] This question concerns Bayesian confidence intervals: \begin{itemize} \item [(i)] Assuming the informative prior used in assignment 1 (not the ignorance prior), which of the following represent valid 95\% posterior Bayesian confidence intervals for $\theta$? Explain your answer. \textcolor{gray}{\bf [3 marks]} \[ (0.2463, 0.5331) \qquad (0.1342, 0.4359) \qquad (0.1102, 0.5021) \qquad (0.1422, 0.4461) \] \item [(ii)] Which interval listed in part (a) is a 95\% posterior highest density interval for $\theta$? Explain your answer. \textcolor{gray}{\bf [2 marks]} \item [(iii)] Find the 95\% prior highest density interval for $\theta$, and compare it to that obtained from the posterior in part (ii). \textcolor{gray}{\bf [2 marks]} \end{itemize} \item [(b)] This question is about the Bayesian posterior predictive distribution: \begin{itemize} \item [(i)] For the next ten major dust storms to hit Greece and Cyprus this year, find the predictive probability distribution for the number $Y$ that will result in an environmental health warning. \textcolor{gray}{\bf [9 marks]} \textcolor{gray}{ Hint: The random variable $X$ follows a Beta-binomial $BetaBin(n, a, b)$ distribution ($n$ positive integer, $a>0, b>0$) if it has probability function \[ f(x|n, a, b) = \left(\begin{array}{c} n\\ x \end{array} \right) \frac{B(x+a, b+n-x)}{B(a,b)}, \quad x=0, 1, \ldots, n, \] where $B(a, b)$ is the beta function and \[ B(a, b) = \frac{(a-1)!(b-1)!}{(a+b-1)!}, \quad \text{for $a$ and $b$ positive integers.} \]} \item [(ii)] Use your answer to part (i) to find a 95\% prediction set for the number of dust storms (of the next ten observed) to trigger an environmental health warning. \textcolor{gray}{\bf [4 marks]} \end{itemize} \end{itemize} \subsection*{Question 2} The negative binomial distribution is used to model scenarios in which we observe the number of independent (success-fail) trials needed before we see a successful trial; the trials must have the same success probability. Suppose \exxvec are a random sample from a negative binomial $NegBin(k,\theta)$ distribution, where $k$ is known. \begin{itemize} \item[(i)] Verify that the congugate prior distribution is a Beta distribution. \item[(ii)] Determine the choice of parameters $g$ and $h$ for the $Beta(g,h)$ distribution that give it maximal variance.\\ Hint: reparameterise the distribution in terms of its mean $m=g/(g+h)$ and $s=g+h$; determine the choice of $m$ and $s$ that maximises the variance of the distribution, and hence the choice of $g$ and $h$. \end{itemize} Determine the posterior distribution for $\theta$ assuming \begin{itemize} \item[(iii)] vague prior knowledge; \item[(iv)] a very large sample. \end{itemize} \subsection*{Question 3} The Rayleigh distribution is often used to measure variability in magnetic resonance imaging (MRI). Suppose \exxvec are a random sample from a Rayleigh $R(\theta)$ distribution. Determine the posterior distribution for $\theta$ assuming \begin{enumerate}[(i)] \item vague prior knowledge; \item a very large sample. \end{enumerate} \subsection*{Question 4} There were 4 ``significant'' flooding events in the UK in 2022. The posterior predictive distribution for $Y$, the number of significant flooding events in 2023, is shown below: \vspace{-0.5cm} \begin{small} \begin{center} \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|}\hline $y$ &0 &1 &2 &3 &4 &5 &6 &7 &8 &9 &10 \\ \hline $f(y|x=4)$ &0.00 &0.01&0.04&0.10 &0.13 &0.45 &0.17 &0.05 &0.03 &0.01 &0.01 \\ \hline \end{tabular} \end{center} \end{small} \begin{itemize} \item [(a)] Find the most plausible 90\% prediction set for $Y$. \item [(b)] Find $E(y|x=4)$ and $SD(y|x=4)$. \end{itemize} \newpage \section*{Solutions} \subsection*{Question 1} \begin{itemize} \item [(a)] Bayesian confidence intervals: \begin{itemize} \item [(i)] Recall from question 1 in assignment 1 that the posterior for $\theta$ is $\theta|\boldsymbol{x} \sim Be(9,23)$. A valid 95\% Bayesian confidence interval $(a, b)$ for $\theta|\boldsymbol{x}$ would be such that $Pr(a<\theta pbeta(0.5331, 9, 23) - dbeta(0.2463, 9, 23) [1] 0.6501951 > pbeta(0.4359, 9, 23) - dbeta(0.1342, 9, 23) [1] 0.9499906 > pbeta(0.5021, 9, 23) - dbeta(0.1102, 9, 23) [1] 0.9899973 > pbeta(0.4461, 9, 23) - dbeta(0.1422, 9, 23) [1] 0.9500078 \end{verbatim} Thus, both $(0.1342, 0.4359)$ and $(0.1422, 0.4461)$ are valid 95\% Bayesian confidence intervals. \textcolor{red}{[1 mark for correct identification, 2 marks for reasoning]} \item [(ii)] We know from the course that a highest density interval also has equal density ($\gamma$) at both $a$ and $b$. Again, we can check in R (or we could do so by hand!). For the first interval: \begin{verbatim} > dbeta(0.1342, 9, 23) [1] 0.8015702 > dbeta(0.4359, 9, 23) [1] 0.8010577 \end{verbatim} And for the second interval: \begin{verbatim} > dbeta(0.1422, 9, 23) [1] 1.038543 > dbeta(0.4461, 9, 23) [1] 0.6451823 \end{verbatim} Thus, only the first interval $(0.1342, 0.4359)$ satisfies $\pi(\theta=a|\boldsymbol{x})=\pi(\theta=b|\boldsymbol{x})=\gamma$, and so $(0.1342, 0.4359)$ is the highest density interval. \textcolor{red}{[2 marks]} We could check this using the \texttt{nclbayes} package: \begin{verbatim} > hdiBeta(0.95, 9, 23) [1] 0.1341826 0.4358943 \end{verbatim} \item [(iii)] We could use the \texttt{nclbayes} package: \begin{verbatim} > hdiBeta(0.95, 5, 17) [1] 0.06921858 0.39949020 \end{verbatim} Thus, the prior HDI is $(0.069, 0.399)$ \textcolor{red}{[1 mark]}-- similar to that for the posterior (it overlaps!), but shifted to the left -- suggesting we are now slightly more inclined to believe that a dust storm will trigger a health warning. \textcolor{red}{[1 mark]} \end{itemize} %\newpage \item [(b)] Predictive distribution: \begin{itemize} \item [(i)] Suppose the number of dust storms triggering a warning in the future is $Y|\theta\sim Bin(10,\theta)$. The predictive probability function of $Y$ is, for $y=0,1,2,\ldots,10,$ \begin{eqnarray*} f(y|x=4) &=&\int_\Theta f(y|\theta)\,\pi(\theta|x=4)\,d\theta \\ & & \\ &=&\int_0^1 \binom{10}{y}\theta^y(1-\theta)^{10-y} \times \frac{\theta^{8}(1-\theta)^{22}}{B(9,23)}\,d\theta \\ & & \\ &=&\binom{10}{y}\frac{1}{B(9,23)}\int_0^1 \theta^{y+8}(1-\theta)^{32-y}\,d\theta \quad \text{\textcolor{red}{[2 marks]}}\\ & & \\ &=& \binom{10}{y}\frac{\color{blue}B(y+9, 33-y)}{B(9,23)} \int_0^1 \frac{\theta^{(y+9)-1}(1-\theta)^{(33-y)-1}}{{\color{blue}B(y+9, 33-y)}}\,d\theta \quad \text{\textcolor{red}{[4 marks]} }\\ &=&\binom{10}{y}\frac{B(y+9,33-y)}{B(9,23)}. \quad \text{\textcolor{red}{[1 mark]}} \end{eqnarray*} That is, $Y|x=4\sim BetaBin(10,9,23)$. \textcolor{red}{[2 mark]} \item [(ii)] Probably best to write an R function? \begin{verbatim} > betabin = function(y, n, a, b){ numerator = beta(y + a, n + b -y) denominator = beta(a, b) answer = choose(10,y)*numerator/denominator return(answer)} > cbind(0:10,round(betabin(0:10, 10, 9, 23),3)) [,1] [,2] [1,] 0 0.058 [2,] 1 0.162 [3,] 2 0.235 [4,] 3 0.230 [5,] 4 0.166 [6,] 5 0.093 [7,] 6 0.040 [8,] 7 0.013 [9,] 8 0.003 [10,] 9 0.000 [11,] 10 0.000 \end{verbatim} The above is our predictive distribution, that is, $f(y|x=4)$. For the prediction set, we need to include values for $y$ with greatest predictive probability, such that the sum of their probabilities is close to 0.95. We include, in order (and with cumulative probabilities in parentheses): $[(2(0.235), 3(0.465), 4(0.631), 1(0.793), 5(0.886), 0(0.944)]$, i.e. $[0, 5]$. \textcolor{red}{[3 marks for method, 1 mark for final answer of $[0, 5]$.} \end{itemize} \end{itemize} \end{document} Therefore, \begin{eqnarray*} f(Y|x=4) &=&\binom{10}{y}\frac{B(y+9,33-y)}{B(9,23)} \\ & & \\ &=& \binom{10}{y} \frac{\frac{(y+9-1)!(33-y-1)!}{(y+9+33-y-1)!}}{\frac{(2-1)!(23-1)!}{(2+23-1)!}} \pause \\ & & \\ &=&\binom{5}{y} \times \frac{(y+1)!(27-y)!\times 24!}{22! \times 29!}. \end{eqnarray*} \pause So, for example, \[ P(Y=0|x=1) = \binom{5}{0} \times \frac{1! \times 27! \times 24!}{22! \times 29!} = 0.680 \quad \text{(See Table 4.4)} \] } } \section*{Distributions for data} \subsection*{Binomial distribution} If $X|\theta\sim Bin(k,\theta)$ then it has probability function \[ f(x|\theta)=\begin{pmatrix}k\\x\end{pmatrix}\theta^x(1-\theta)^{k-x} \,,\quad x=0,1,\ldots,k, \] where $k$ is a positive integer and $0<\theta<1$. Also, $E(X)=k\theta$ and $Var(X)=k\theta(1-\theta)$. \subsection*{Exponential distribution} If $X|\lambda\sim Exp(\lambda)$ then it has density \[ f(x|\lambda)=\lambda e^{-\lambda x}\,,\quad x>0 \] where $\lambda>0$. Also, $E(X)=1/\lambda$ and $Var(X)=1/\lambda^2$. \subsection*{Gamma distribution} If $X|\alpha,\lambda\sim Ga(\alpha,\lambda)$ then it has density \[ f(x|\alpha,\lambda)=\frac{\lambda^\alpha x^{\alpha-1}e^{-\lambda x}}{\Gamma(\alpha)}\,, \quad x>0, \] where $\alpha>0$ and $\lambda>0$. Also, $E(X)=\alpha/\lambda$ and $Var(X)=\alpha/\lambda^2$. \subsection*{Inverse Gaussian distribution} If $X|\mu,\lambda\sim IG(\mu,\lambda)$ then it has density \[ f(x|\mu,\lambda)=\frac{\sqrt\lambda}{\sqrt{2\pi x^3}} \exp\left\{-\frac{\lambda(x-\mu)^2}{2\mu^2 x}\right\}\,, \quad x>0, \] where $\mu>0$ and $\lambda>0$. Also, $E(X)=\mu$ and $Var(X)=\mu^3/\lambda$. \subsection*{Laplace distribution} If $X|\mu,\sigma\sim La(\mu,\sigma^2)$ then it has density \[ f(x|\mu,\sigma)=\frac{1}{\sqrt{2}\,\sigma} \exp\left\{-\frac{\sqrt{2}|x-\mu|}{\sigma}\right\}\,, \quad x>0, \] where $\mu\in\mathbb{R}$ and $\sigma>0$. Also, $E(X)=\mu$ and $Var(X)=\sigma^2$. \subsection*{Log-normal distribution} If $X|\mu,\sigma\sim LN(\mu,\sigma^2)$ then it has density \[ f(x|\mu,\sigma)= \frac{1}{x\sigma\sqrt{2\pi}} \exp\left\{-\frac{1}{2\sigma^2}(\log x-\mu)^2\right\},\quad x>0 \] where $\mu\in\mathbb{R}$ and $\sigma>0$. Also, $E(X)=e^{\mu+\sigma^2/2}$, $Var(X)=(e^{\sigma^2}-1)e^{2\mu+\sigma^2}$. Further $E(\log X)=\mu$ and $Var(\log X)=\sigma^2$. \subsection*{Negative binomial distribution} If $X|\theta\sim NegBin(k,\theta)$ then it has probability function \[ f(x|\theta)=\binom{x-1}{k-1}\theta^k(1-\theta)^{x-k},\quad x=k,k+1,\ldots~, \] where $k$ is a positive integer and $0<\theta<1$. Also $E(X)=k/\theta$ and $Var(X)=k(1-\theta)/\theta^2$. \subsection*{Normal distribution} If $X|\mu,\tau\sim N(\mu,1/\tau)$ then it has density \[ f(x|\mu,\tau)=\left(\frac{\tau}{2\pi}\right)^{1/2} \exp\left\{-\frac{\tau}{2}(x-\mu)^2\right\},\quad x\in\mathbb{R} \] where $\mu\in\mathbb{R}$ and $\tau>0$. Also, $E(X)=\mu$ and $Var(X)=1/\tau$. The distribution has the following quantiles \begin{center} \begin{tabular}{c|ccccc} $ x $ & 1.2816 & 1.6449 & 1.9600 & 2.3263 & 2.5758 \\ \hline $ \Pr(X\theta \] where $\alpha>0$ and $\theta>0$. Also $E(X)=\alpha\theta/(\alpha-1)$, $\alpha>1$ and $Var(X)=\alpha\theta^2/\{(\alpha-1)^2(\alpha-2)\}$, $\alpha>2$. \subsection*{Poisson distribution} If $X|\theta\sim Po(\theta)$ then it has probability function \[ f(x|\theta)=\displaystyle \frac{\theta^xe^{-\theta}}{x!}\,, \quad x=0,1,\ldots, \] where $\theta>0$. Also, $E(X)=\theta$ and $Var(X)=\theta$. \subsection*{Rayleigh distribution} If $X|\theta\sim R(\theta)$ then it has density \[ f(x|\theta)=2x\theta e^{-\theta x^2},\quad x>0, \] where $\theta>0$. Also $E(X)=\sqrt{\pi/(4\theta)}$ and $Var(X)=(4-\pi)/(4\theta)$. \subsection*{Uniform distribution} If $X|\phi,\theta\sim U(\phi,\theta)$ then it has density \[ f(x|\phi,\theta)=\frac{1}{\theta-\phi}\,,\quad \phi0$. Also $I_0(\cdot)$ is the modified Bessel function of order zero. \subsection*{Weibull distribution} If $X|\beta,\lambda\sim Wei(\beta,\lambda)$ then it has density \[ f(x|\beta,\lambda)=\beta\lambda^\beta x^{\beta-1}e^{-\lambda^\beta x^\beta}\,, \quad x>0, \] where $\beta>0$ and $\lambda>0$. Also, $E(X)=\Gamma(1+1/\beta)/\lambda$ and $Var(X)=\Gamma(1+2/\beta)/\lambda^2-E(X)^2$. \newpage \section*{Distributions for prior beliefs} \subsection*{Beta distribution} If $\theta\sim Beta(g,h)$ then it has density \[ \pi(\theta)=\frac{\theta^{g-1}(1-\theta)^{h-1}}{B(g,h)}\,, \quad 0<\theta<1, \] where $g>0$ and $h>0$. Also, $E(\theta)=g/(g+h)$ and $Var(\theta)=gh/\{(g+h)^2(g+h+1)\}$. \subsection*{Exponential distribution} If $\theta\sim Exp(h)$ then it has density \[ \pi(\theta)=h e^{-h\theta}\,,\quad \theta>0, \] where $h>0$. Also, $E(\theta)=1/h$ and $Var(\theta)=1/h^2$. \subsection*{Gamma distribution} If $\theta\sim Ga(g,h)$ then it has density \[ \pi(\theta)=\frac{h^g \theta^{g-1}e^{-h\theta}}{\Gamma(g)}\,, \quad \theta>0, \] where $g>0$ and $h>0$. Also, $E(\theta)=g/h$ and $Var(\theta)=g/h^2$. \subsection*{Generalised $\boldsymbol{t}$ distribution} \label{page:gt} If $\mu\sim t_a(b,c)$ then it has density \[ \pi(\mu)=\frac{\Gamma\left(\frac{a+1}{2}\right)} {\sqrt{ac\pi}\,\Gamma\left(\frac{a}{2}\right)} \left\{1+\frac{(\mu-b)^2}{ac}\right\}^{-\frac{a+1}{2}},\quad \mu\in\mathbb{R}, \] where $b\in\mathbb{R}$, $a>0$ and $c>0$. Also, $E(\mu)=b$ and $Var(\mu)=ac/(a-2)$ if $a\geq 2$. \subsection*{Inverse Chi distribution} \label{page:invchi} If $\sigma\sim\textit{Inv-Chi}(a,b)$ then it has density \[ \pi(\sigma|a,b)=\frac{2b^a\sigma^{-2a-1}e^{-b/\sigma^2}}{\Gamma(a)},\quad \sigma>0, \] where $a>0$, $b>0$ and $\Gamma(a)$ is the gamma function. Also $E(\sigma)=\sqrt{b}\,\Gamma(a-1/2)/\Gamma(a)$ and $Var(\sigma)=b/(a-1)-E(\sigma)^2$ if $a>1$. The name of the distribution comes from the fact that $1/\sigma^2\sim Ga(a,b)\equiv \chi^2_{2a}/(2b)$. \subsection*{Log-normal distribution} If $\theta\sim LN(b,c^2)$ then it has density \[ \pi(\theta)=\frac{1}{\sqrt{2\pi}\,c\,\theta} \exp\left\{-\frac{1}{2c^2}(\log\theta-b)^2\right\},\quad \theta>0 \] where $b\in\mathbb{R}$ and $c>0$. Also, $E(\theta)=e^{b+c^2/2}$, $Var(\theta)=(e^{c^2}-1)e^{2b+c^2}$. Further $\log\theta\sim N(b,c^2)$ and so $E(\log\theta)=b$ and $Var(\log\theta)=c^2$. \subsection*{Normal distribution} If $\mu\sim N(b,1/d)$ then it has density \[ \pi(\mu)=\left(\frac{d}{2\pi}\right)^{1/2} \exp\left\{-\frac{d}{2}(\mu-b)^2\right\},\quad \mu\in\mathbb{R}, \] where $b\in\mathbb{R}$ and $c>0$. Also, $E(\mu)=b$ and $Var(\mu)=1/d$. \subsection*{Normal-gamma distribution} If $\mutau\sim NGa(b,c,g,h)$ then it has density \[ \pi(\mu,\tau)\propto\tau^{g-\frac{1}{2}}\exp\left\{-\frac{\tau}{2} \left[c(\mu-b)^2+2h\right]\right\},\quad \mu\in\mathbb{R}, ~\tau>0 \] where $b\in\mathbb{R}$ and $c,g,h>0$. Also, $\mu|\tau\sim N\left(b,\dfrac{1}{c\tau}\right)$, $\tau\sim Ga(g,h)$ and has marginal distribution $\mu\sim t_{2g}\left(b,\dfrac{h}{gc}\right)$. \subsection*{Uniform distribution} If $\theta\sim U(a,b)$ then it has density \[ \pi(\theta)=\frac{1}{b-a}\,,\quad a<\theta