MATH454/554: Project II

TO BE HANDED IN BY MONDAY 11/12/2017 (WEEK 10), 10:00.

This project will contribute 17% towards the final mark.

Submission: Upload the pdf of your answer and your R code file to the Moodle site. Your R code should be as .r or .txt file so that it can be copied and pasted to run. Submit also a printed copy of your answers (no need for a printed copy of the R code), together with a plagiarism cover sheet, to the MSc submissions pigeon hole. Please write your student ID on your answers, not your name.

This project looks at modelling the number of movements of a fetal lamb over 240 consecutive periods, each of 5 seconds. The data which is given in lamb.txt comes from Leroux and Puterman (1992) and has also been analysed in Fearnhead (2005).

Let $𝐱=(x_1,x_2,\mathrm{\dots },x_{240})$ denote the 240 fetal lamb movements. It is assumed that the data arise from the following Markov-dependent mixture model defined in Fearnhead (2005), Section 3.3. There are two underlying (unobserved) states, which we will term states 1 and 2. Let $y_t$ denote the state at time $t$. Then it is assumed that

Furthermore, it is assumed that the unobserved states $𝐲=(y_1,y_2,\mathrm{\dots },y_{240})$ follow a Markovian structure with $y_1=1$ and for $t>1$,

$P(y_t=1|y_{t-1}=1)$ $=$ $p_1$

$P(y_t=2|y_{t-1}=1)$ $=$ $1-p_1$

$P(y_t=1|y_{t-1}=2)$ $=$ $1-p_2$

$P(y_t=2|y_{t-1}=2)$ $=$ $p_2.$

Thus the process $\{y_t\}$ is Markov chain with transition matrix

We are interested in constructing a Gibbs sampler to obtain samples from $\pi (𝜽|𝐱)$, where $𝜽=({\lambda }_1,{\lambda }_2,p_1,p_2)$. We will use data augmentation of $𝐲$ to achieve this. (Note that $y_1=1$ is known and therefore fixed.)

Assume independent priors for the four parameters with a $\mathrm{Gamma}(1,1)$ prior for ${\lambda }_1$ and ${\lambda }_2$ and a $\mathrm{Beta}(1,1)$ prior (uniform prior) for $p_1$ and $p_2$.

1. Write down $\pi (𝐱,𝐲|𝜽)$, the likelihood of the observed data (number of movements of the fetal lamb) and the augmented data (mixture component). [2]

2. For , show that [1]

   $$\pi (y_t|𝐱,{𝐲}_{-t},𝜽)\propto \frac{{\lambda }_{y_t}^{x_t}}{x_t!}\exp (-{\lambda }_{y_t})p_{y_{t-1}}^{I(y_t=y_{t-1})}{(1-p_{y_{t-1}})}^{I(y_t\ne y_{t-1})}{p}_{y_t}^{I(y_{t+1}=y_t)}{(1-p_{y_t})}^{I(y_{t+1}\ne y_t)}.$$

   Hence, show that [2]

   $$P(y_t=1|𝐱,{𝐲}_{-t},𝜽)=\frac{Q(1;y_{t-1},y_{t+1},x_t,𝜽)}{Q(1;y_{t-1},y_{t+1},x_t,𝜽)+Q(2;y_{t-1},y_{t+1},x_t,𝜽)},$$

   where for $s=1,2$,

   $$Q(s;y_{t-1},y_{t+1},x_t,𝜽)={\lambda }_s^{x_t}\exp (-{\lambda }_s)p_{y_{t-1}}^{I(s=y_{t-1})}{(1-p_{y_{t-1}})}^{I(s\ne y_{t-1})}{p}_s^{I(y_{t+1}=s)}{(1-p_s)}^{I(y_{t+1}\ne s)}.$$

   Similarly,

   $P(y_{240}=1|𝐱,{𝐲}_{-240},𝜽)$

   $=$ ${\displaystyle \frac{{\lambda }_1^{x_{240}}\exp (-{\lambda }_1)p_{y_{239}}^{I(1=y_{239})}{(1-p_{y_{239}})}^{I(1\ne y_{239})}}{{\lambda }_1^{x_{240}}\exp (-{\lambda }_1)p_{y_{239}}^{I(1=y_{239})}{(1-p_{y_{239}})}^{I(1\ne y_{239})}+{\lambda }_2^{x_{240}}\exp (-{\lambda }_2)p_{y_{239}}^{I(2=y_{239})}{(1-p_{y_{239}})}^{I(2\ne y_{239})}}}.$

3. For $j=1,2$, compute the conditional distribution of ${\lambda }_j$ given $𝐱,𝐲$ and ${𝜽}_{-{\lambda }_j}$. [2]

4. For $j,k=1,2$, let $M_{jk}={\sum }_{t=2}^{240}I(y_{t-1}=j)I(y_t=k)$. Show that for $j=1,2$, that the conditional distribution of $p_j$ given $𝐱,𝐲$ and ${𝜽}_{-p_j}$ satisfies [2]

   $$p_j|𝐱,𝐲,{𝜽}_{-p_j}\sim \mathrm{Beta}(M_{jj}+1,M_{j(3-j)}+1).$$

5. Write an R routine to implement the Gibbs sampler. [5]

6. Run the R routine to obtain a sample of size 51000 from the posterior and discard the first 1000 iterations as burn-in. Compute the posterior means and standard deviations of the parameters. [3]