Project 2

Submit here no later 23:59 on Sunday March 16.

Oral presentations at 16:15 on Thursday March 13. Submit you presentation here. The plan for the presentations (which group presents what) is here.

Project supervision on Tuesdays 10:15-12:00 in weeks 10 and 11.

Further guidelines are given here (same as for project 1).

Reconsider problem 3a from project 1, that is, how to sample from a bivariate standard normal with zero correlation truncated to a disc with radius \(\rho\) centered at \(\delta,0\). The aim here is to construct a simple Gibbs sampling algorithm for this problem.

a) Write down an expression for the joint target density function up to a normalising constant in Cartesian coordinates. Optionally, you may want to write down the joint density after transforming to polar coordinates. Working either with Cartesian or polar coordinates, derive expressions for the full conditional densities (including their normalising constants) of each coordinate (that is, either \(f_{X|Y}(x|y)\) and \(f_{Y|X}(y|x)\) or \(f_{R|\Theta}(r|\theta)\) and \(f_{\Theta|R}(\theta|r)\)) and derive inversion methods for simulating from each full conditional. Hint: Recall that the full conditionals are always proportional to the joint density. To understand the geometry of the problem you will want to make a drawing. You may also find either the law of cosines or this equation for the \(y\)-coordinate of a circle as function of \(x\) useful.

b) Implement an R-function that simulates \(n\) realisations from the target density using Gibbs sampling. The function should return the realisations in Cartesian coordinates as a \(n \times 2\) matrix of class mcmc as defined by coda R-package. Examine trace plots (obtain by applying plot to the returned mcmc object) of each variable and briefly comment on what you see.

c) Use your Gibbs sampler to generate \(10^4\) samples for \(\delta=2.5\) and \(\rho=0.5\). Using the effectiveSize function in the coda package, compute an estimate the effective sample size (use the estimate for the coordinate with smallest size). Using the microbenchmark package, measure CPU-time consumed by your Gibbs sampler. Also run the expression

 microbenchmark::microbenchmark({n <- 100; x <- solve(matrix(runif(n^2),n,n))})

Report your measurements here.

Suppose that we observe a renewal process over a closed time interval \([0,t]\) and that \(n\) failures are observed occuring at times \(y_1,y_2,\dots,y_n\) where \(0<y_1<y_2<\dots<y_n<t\). We will assume that the inter-arrival times are independent and Gamma distributed with shape parameter \(\alpha\) and rate parameter \(\lambda\), that is, the density of \(y_i\) conditional on \(y_{i-1}\) is given by \[ f(y_i|y_{i-1})=\frac{\lambda^\alpha}{\Gamma(\alpha)}(y_i-y_{i-1})^{\alpha-1}e^{-\lambda(y_i-y_{i-1})}I(y_i>y_{i-1}) \] for \(i=2,3,\dots,n\). Since we begin observing the process at an arbitrary time (without loss of generality set to zero), it can be shown that the distribution of \(y_1\) is not Gamma but instead follows the residual time density \[ \phi(y_1)=\frac{1-F(y_1)}\mu I(y_1\ge 0) \] where \(F\) is the cdf of the interarrival times and \(\mu\) is their expected value (thus given by \(\mu=\alpha/\lambda\) for the above model).

We will represent our prior beliefs about \(\alpha\) and \(\lambda\) by the prior \[ \pi(\lambda,\alpha)\propto \frac1\lambda e^{-\alpha}I(\alpha>0)I(\lambda>0), \] that is, the shape is exponentially distributed with a prior expected value of 1 and the rate parameter has an improper scale prior (see e.g. Berger (1985), p. 83 for a rationale behind this prior).

The observed data are given by

y <- c(5.14, 10.15, 18.4, 31.72, 35.4, 36.5, 37.95, 40.82, 64.84, 
96.8)
t <- 100

a) Given the model parameters, what is the probability of observing no further arrivals on \([0,t]\) after \(y_n\)? Write down an expression for the joint posterior density of \(\alpha,\lambda\) up to a normalising constant. Hint: Use the probability chain rule. As a sanity check, how does your expression simplify in the special case of \(\alpha=1\)?

To sample from the posterior density of \(\alpha,\lambda\) we will use Gibbs sampling. To this end, we will augment the data with two additional latent variables \(y_0<0\) and \(y_{n+1}>t\) representing the unknown arrival times just before and after the beginning and end of the observation time interval \([0,t]\).

b) What are the full conditional densities of \(y_0\), \(y_{n+1}\), that is, conditional on the observed data and \(\alpha\) and \(\lambda\)? Hint: A renewal process is time-reversible. How can we simulate from these conditional densities using the inversion method?

c) Derive the joint full conditional of \((\alpha,\lambda)\) up to a normalising constant. Hint: When we condition also on \(y_0\) and \(y_{n+1}\) we essentially have \(n+2\) observations on the closed interval \([y_0,y_{n+1}]\) which leads to an expression for the joint posterior similar to the expression in point a) but with some terms simplifying. Show that full conditional of \(\lambda\) only now is another Gamma distribution. Also write down the full conditional of \(\alpha\) only up to a normalising constant.

d) Implement a Gibbs sampler (as a function gibbs in R) to obtain samples from the posterior distribution of \(\alpha,\lambda,y_0,y_{n+1}\) updating \(\lambda\), \(\alpha\), \(y_0\) and \(y_{n+1}\) in separate Gibbs steps. You will need to use Metropolis-within-Gibbs for \(\alpha\), e.g. a Gaussian random-walk proposal. Your function should return the realised Markov-Chain as an object of class mcmc as defined by the coda R-package. Examine trace plots for each parameter. Also make a scatter plots of the samples for \(\lambda\) and \(\alpha\). Tune the step size of the random-walk proposal for \(\alpha\) to improve convergence. Brifely discuss why the convergence of the Markov-chain is rather slow.

e) To improve convergence, \((\alpha,\lambda)\) can be updated in a single Metropolis-within-Gibbs step using a block proposal of the form \(Q(\alpha',\lambda'|\alpha,\lambda,\mathbf{y})=Q(\alpha'|\alpha)\pi(\lambda'|\alpha',\mathbf{y})\) where \(Q(\alpha'|\alpha)\) is the random-walk proposal in point d). Briefly discuss the intuition behind why this might be more efficient. Derive a fully simplified expression for the log of the resulting acceptance probability and create a new function gibbsblock implementing this strategy. Tune the variance of random-walk proposal to improve the rate of convergence of the resulting chain. Compare the behaviour of the resulting Markov-chain (including the autocorrelation and effective sample size of the resulting chain) with the algorithm in point d).

Time series measurements are often observed with data irregularities, e.g., observations with a detection limit. For instance, a monitoring device usually has a detection limit and it records the limit value when the true value exceeds/precedes the detection limit. In this problem you will analyse a time series of cloud ceiling heights taken from a paper by Park et.al. (2007), collected by the National Center for Atmospheric Research (NCAR), observed hourly in San Francisco and recorded during the month of March 1989, consisting of n = 716 observations. Cloud ceiling height is one of the major factors contributing to weather-related accidents and one of the major causes of flight delays. The observations have a detection limit at about 4.79 (12,000 feet in original scale) and hence we consider the data as a right-censored time series. Due to this upper detection limit of the record- ing device, many observations are censored: the censoring rate is 41%.

Load the dataframe CloudCeiling into your R-workspace and plot the data using the command

 load(url("https://www.math.ntnu.no/emner/TMA4300/2025v/CloudCeiling.RData"))     
 plot(CloudCeiling$y, xlab="Time", ylab="log(Cloud ceiling height)")
 

The variable y contains the logarithm of hourly cloud ceiling heights (in units of 100 feets on the original scale), and cc is a right censoring indicator (1 if the observation is right-censored, 2 if the observation is completely missing, and 0 otherwise). Missing observations (of which there are 3 have been imputed preliminary values by interpolation).

Using Bayesian inference, our aim is to estimate the parameters of a first order auto-regressive model of the form \[ y_t - \mu = \phi(y_{t-1} - \mu) + \epsilon_t \tag{1} \] where \(\epsilon_t\) is zero-mean iid Gaussian white noise with variance \(\tau^{-1}\) as well as the posterior density of the censored observations. As joint prior we will use the improper density \[ \pi(\mu, \phi, \tau) \propto I(-1<\phi<1)\frac1{\tau}. \] Thus we are assuming a priori that the process is stationary.

a) Show that stationarity implies that \(E(y_1|\mu,\color{red}{\phi},\tau)=\mu\) and that \(\operatorname{Var}(y_1|\mu,\color{red}{\phi},\tau)=\frac{\tau^{-1}}{1-\phi^2}\). Corrected on March 11.

At first we will ignore the fact that we have right-censored (and some missing) data and treat the censoring times of 4.78 (and interpolated values) as if they are observed cloud ceiling heights.

b) Reparameterizing the model by writing it on the form \[ y_t = \beta_0 + \beta_1 y_{t-1} + \epsilon_t, \tag{2} \] show that the joint prior on \(\beta_0, \beta_1, \tau\) becomes \[ \pi(\beta_0, \beta_1, \tau) = I(-1<\beta_1<1)\frac1{1-\beta_1}\frac1{\tau}. \] Letting \(\boldsymbol\beta=(\beta_0, \beta_1)^T\) and \(\mathbf{X} = \begin{bmatrix} \mathbf{1} & \mathbf{y}_{1:(n-1)}\end{bmatrix} \) show that the joint posterior of the model parameters can be written as \[ \pi(\boldsymbol\beta,\tau|\mathbf{y}) \propto \underbrace{\tau^{\frac n2-1} \exp\Big(-\frac\tau 2(\mathbf{y}_{2:n}-\mathbf{X}\boldsymbol\beta)^T (\mathbf{y}_{2:n}-\mathbf{X}\boldsymbol\beta)\Big)}_{(*)}\frac{\sqrt{1-\beta_1^2}I(-1<\beta_1<1)}{1-\beta_1}\exp\Big(-\frac{\tau}2(1-\beta_1^2)(y_1-\frac{\beta_0}{1-\beta_1})^2\Big) \]

c) The part denoted (*) has the same form (up to a normalising constant) as the posterior density of an ordinary linear regression model with non-informative prior. If we use (*) as a joint proposal density \(Q(\boldsymbol\beta',\tau'|\boldsymbol\beta,\tau,\mathbf{y})\) in a Metropolis-Hastings (independence sampler) algorithm, derive a fully simplified expression for the log of the resulting acceptance probability \(\alpha(\boldsymbol\beta',\tau'|\boldsymbol\beta,\tau,\mathbf{y})\).

d) Show that (*), up to a constant of proportionality, factors into \[ |\mathbf{Q}|^{\frac12}\tau^{\frac p2}\exp\Big(-\frac\tau2(\boldsymbol\beta-\hat{\boldsymbol\beta})^T \mathbf{Q} (\boldsymbol\beta-\hat{\boldsymbol\beta})\Big)\tau^{\frac{n-p}2-1} \exp\Big(-\frac{\tau}2 \hat{\boldsymbol\epsilon}^T\hat{\boldsymbol\epsilon} \Big) \] where \(p=2\) in our case, that is a multivariate normal density for \(\boldsymbol\beta|\tau\) with precision matrix \(\tau\mathbf{Q}=\tau\mathbf{X}^T \mathbf{X}\) and mean vector \(\hat{\boldsymbol\beta}=(\mathbf{X}^T \mathbf{X})^{-1}\mathbf{X}^T \mathbf{y}_{2:n}\) and a Gamma distribution for the marginal density of \(\tau\) with shape parameter \(\frac{n-p}2\) and rate parameter \(\frac12 \hat{\boldsymbol\epsilon}^T\hat{\boldsymbol\epsilon}\) where \(\hat{\boldsymbol\epsilon}=\mathbf{y}_{2:n}-\mathbf{X}\hat{\boldsymbol\beta}\). Explain how this can be exploited to simulate from the joint proposal in point c.

e) Implement a Metropolis-Hastings algorithm (an independence sampler) as a function in R for simulating from \(\pi(\boldsymbol\beta,\tau|\mathbf{y})\). You may want to use the rmvnorm function in the mvtnorm R-package. Compute the overall acceptance rate and examine the convergence of the chain. Why is (*) a good or a bad choice of proposal?

Having implemented a method treating all the data as non-censored, we now turn to how to account for the right censoring of some of the observations. For observations that are right censored or missing, we will now instead let \(y_i\) represent latent unobserved cloud ceiling heights for these observations.

f) Let \(\mathbf{y}_{-j}\) denote the subvector obtained by removing \(y_j\) from \(\mathbf{y}\) and let \(\mathbf{c}=(c_1,c_2,\dots,c_n)\) denote the censoring indicators. In terms of parameterisation (1), show that \[ y_j|\mathbf{y}_{-j},\boldsymbol\beta,\tau \sim N\left(\mu + \frac\phi{1+\phi^2}(y_{j-1}+y_{j+1}-2\mu),\frac{\tau^{-1}}{1+\phi^2}\right). \] What is the density of \(y_j|\mathbf{y}_{-j},\mathbf{c},\boldsymbol\beta,\tau\) in cases where the censoring indicator \(c_j=1\)?

Hint: It follows that \[ \begin{bmatrix}y_{j-1} \\ y_j \\ y_{j+1}\end{bmatrix} \sim N\left( \mu\mathbf{1}, \tau^{-1} \begin{bmatrix} 1 & -\phi & 0 \\ -\phi & 1+\phi^2 & -\phi \\ 0 & -\phi & 1 \end{bmatrix}^{-1} \right). \] Also see lecture notes from part 1.

g) Implement a Gibbs sampler updating \(\boldsymbol\beta,\tau\) from its full conditional using the method in point e) and updating each censored or missing \(y_j\) in separate Gibbs steps by simulating from the density in point f). Return the chain (including realisations of censored \(y_j\), as an object of class mcmc. Examine the rate of convergence by studying traceplots, estimating the autocorrelation for the different parameters, and by computing effective sample sizes. Increase the number of iterations accordingly. Plot some realisations of all of \(\mathbf{y}\) and comment on what you see. Discuss how the overall rate of convergence is likely to depend on value of \(\phi\) suggested by the data. How do the parameter estimates change when we account for the right-censoring of the data?

In this problem, using RTMB, we will analyse the seasonal pattern in the probability of rainfall in Tokyo on each day of the year. Load the dataframe rain into your R workspace using the command

 load(file=url("https://www.math.ntnu.no/emner/TMA4300/2025v/rain.rda"))

The variable n.rain (denoted \(y_i\) in the following) contains the number of years out of n.year (denoted \(n_i\)) that it rained on different particular days of the year given by the variable day (denoted \(i\) in the following) taking values from 1 to \(n=366\).

For each day \(i=1,2,\dots,n\) we will assume that \[ y_i|x_i \sim \operatorname{bin}(n_i, p_i) \] where \[ \operatorname{logit}p_i = x_i. \] We will model \(\mathbf x=(x_1,x_2,\dots,x_n)^T\) as an Intrinsic Gaussian Random Markov Field (IGRMF). A random walk is an example of a IGRMF. However, since the seasonal trend in probability of rainfall a priori is expected to be a smooth function, a second order random walk may more reasonable. This amounts to the assumption that the second order differences \[ (1-B)^2 x_t = (1-2B+B)x_t = x_t-2x_{t-1}+x_{t-2} \sim N(0,\kappa^{-1}) \] where \(B\) is the backshift operator. Thus, the parameter \(\kappa\) determines the smoothness of realisations of \(\mathbf{x}\). Under this model, \(\mathbf x\) has a joint improper density given by \[ \pi(\mathbf x|\kappa) \propto \kappa^{\frac{n-2}2}\exp\left\{-\frac\kappa2\sum_{i=3}^n (x_i-2x_{i-1}+x_{i-2})^2\right\} =\kappa^{\frac{n-2}2}\exp\left\{-\frac\kappa2\mathbf{x}^T\mathbf{Q}\mathbf{x}\right\} \] The density is not a proper multivariate Gaussian density since the precision matrix \(\mathbf{Q}\) can be shown to have rank \(n-2\). Thus, the density is flat in two directions. See Rue and Held, 2005, for details. It should be further noted that since the density is improper, the "normalising constant" is arbitrary and there is disagreement about weather \(\frac{n-2}2\) is the correct power of \(\kappa\), see Lavine and Hodges (2012).

a) Under the above assumptions, what is the joint density of the data \(\mathbf{y}\) and the latent variables \(\mathbf{x}\)? Write a R function f computing the negative log of this joint density. The format should be written such that it can be passed as an argument to makeADFun in the RTMB package so you need to study this RTMB introduction and maybe also the R-code from the lecture 2 in week 9, in particular, \(\kappa\) and \(\mathbf{x}\) needs to be passed to f as a list object.

b) What is the definition of the marginal likelihood for the above model and how does it relate to the joint density in point a)? How is the Laplace approximation of the marginal likelihood defined?

c) Run the command obj ← MakeADFun( ) with the necessary arguments including random="x" and maximise the Laplace approximation of the marginal likelihood to compute the MLE of \(\kappa\). Also extract estimates of \(\mathbf{x}\) and associated standard errors using the sdreport and summary functions (again study the RTMB introduction for details). Make a plot of the estimated seasonal trend in both the logit of the probability of rain fall and of the probability itself overlaying the observed data.

d) A shortcoming of the above second order random walk model for \(\mathbf{x}\) is that the estimated trend may have a undesirable discontinuity when \(i\) goes from 365 to 1. This can be solved by assuming that \(\mathbf{x}\) instead follows a cyclic second order random walk which has density \[ \pi(\mathbf x|\kappa) \propto \kappa^{\frac{n-1}2}\exp\left\{-\frac\kappa2\sum_{i=1}^n (x_i-2x_{i-1\operatorname{mod} n)}+x_{i-2\operatorname{mod} n})^2\right\}. \] Note that since we are indexing the observations from 1 to \(n\) rather than 0 to \(n-1\) the Modulo operator must be redefined accordingly such that \(0\operatorname{mod}n\) equals \(n\) instead of 0. This associated circulant precision matrix now instead have rank \(n-1\) suggesting that the correct power of \(\kappa\) is \(\frac{n-1}2\). Change your code such that you instead can fit this model and examine the resulting estimated seasonal trend. How do the maximum (marginal) likelihoods of each model alternative compare?