Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 216 / 3

Write $\lambda ={\sigma }^2$, $v={\sigma }_0^2$, and $a=\sqrt{\rho }$. We keep the model's printed coefficient ${\sigma }^2{Z}_i$ literally. Its latent log odds can be writtenThe independence of the ${\theta }_i$ is conditional on $Z$: they are not marginally independent when the random effects $Z_i$ are correlated.

Holding the other predictors and the random-effect realization fixed, increasing predictor $j$ by one shifts the latent log odds by ${\beta }_j$, so it multiplies the conditional odds of a late departure by $e^{{\beta }_j}$. For a binary snow indicator, this compares the snow and no-snow days with other predictors held fixed. This is a conditional odds ratio, not generally the same as a marginal odds ratio after integrating the random effects.

Randomizing ${\theta }_i$ allows unmeasured daily conditions to change the success probability beyond the systematic linear predictor. In this logistic-normal regression with autoregressive random effects, let $P_i={\mathrm{logit}}^{-1}({\theta }_i)$ and ${\bar{p}}_i=\mathbb{E}[P_i]$. The law of total variance givesso the mixture allows overdispersion beyond a binomial distribution with fixed probability. The correlated part also allows related outcomes on nearby days.

The stationary autoregressive process of order one has $\mathrm{Var}(Z_i)=1$ and $\mathrm{Cov}(Z_i,Z_j)=a^{|i-j|}$. HenceThe ${\sigma }_0^2$ term is unstructured daily variability, whereas ${\sigma }^4$ is the temporally correlated component under the printed parameterization. Larger $\rho$ gives greater persistence: the lag-$h$ correlation of the latent log odds is ${\sigma }^4{\rho }^{h/2}/({\sigma }^4+{\sigma }_0^2)$.

The PDF itself prints both the coefficient ${\sigma }^2{Z}_i$ and the requested combination ${\sigma }^2+{\sigma }_0^2$. They are inconsistent as a marginal variance. If the intended coefficient were $\sigma Z_i$, the variance and off-diagonal covariance would instead be ${\sigma }^2+{\sigma }_0^2$ and ${\sigma }^2{\rho }^{|i-j|/2}$. The prior's argument list also repeats ${\sigma }^2$; its displayed density indicates the intended second variance parameter is ${\sigma }_0^2$.

There is a further substantive issue with that improper prior. After integrating the latent variables, the observed-data likelihood function is positive and continuous at $\lambda =0$ for fixed finite $\beta$, positive $v$, and $\rho \in (0,1)$. On a compact positive-volume set of those parameters it therefore has a positive lower bound for small $\lambda$. The prior density is $1/(\lambda v)$, so integrating over $\lambda$ gives ${\int }_0^{\epsilon }d\lambda /\lambda =\infty$. Thus the stated joint posterior distribution is improper, an instance of improper posterior from a log-uniform random-effect scale prior. The requested fixed-parameter conditional laws below are nevertheless proper; their existence does not remedy the joint impropriety.

Use the printed parameterization $\lambda ={\sigma }^2$, $v={\sigma }_0^2$, and $a=\sqrt{\rho }$. Let $T=365$, $X$ be the design matrix with rows $x_i^T$, and $r=\theta -X\beta$. Conditional on $\theta$ and the parameters, $Y$ contributes no extra dependence on $Z$ by conditional independence.

The stationary autoregressive normal distribution has prior density proportional toIts Gaussian autoregressive conditional precision is the tridiagonal matrix $J$ withCombining this prior with $r\mid Z\sim N(\lambda Z,vI)$ and completing the square givesThe precision matrix $Q$ is positive definite. Compute its lower Cholesky decomposition $Q=L{L}^T$, solve $Qm=h$, draw $\xi \sim N(0,I_T)$, and solve $L^T\eta =\xi$. Then return $Z=m+\eta$. This Gaussian sampling from a precision Cholesky factor is exact because $\mathrm{Cov}(\eta )=L^{-T}{L}^{-1}=Q^{-1}$. Tridiagonal structure gives linear cost in $T$ and preserves the conditional correlation between days. Notice that the solve uses $L^T$, not $L$.

These conditionals are normalized for every fixed positive ${\sigma }^2,{\sigma }_0^2$, despite the joint posterior impropriety explained in part (a). If the intended model uses $\sigma Z_i$ instead, replace the observation coefficient $\lambda$ by $\sigma$ throughout; the information and precision terms then use $\sigma /v$ and ${\sigma }^2/v$.

Given the parameters and $Z$, the coordinates are conditionally independent. Write $m_i=x_i^T\beta +{\sigma }^2{Z}_i$, $v={\sigma }_0^2$, and let $y_i$ denote the observed count. Their full conditional distributions have scalar densitiesEach is proper: the binomial distribution likelihood factor is bounded, and the remaining factor is a proper normal distribution density.

One completely specified exact algorithm is rejection sampling from $N(m_i,v)$. DefineFor $0<y_i<n_i$,for $y_i=0$ or $y_i=n_i$, take $M_i=1$, and for $n_i=0$ take $H_i=M_i=1$. Propose $t\sim N(m_i,v)$ and accept with probability $H_i(t)/M_i$, repeating until acceptance. The accepted density is proportional to the proposal times $H_i$, exactly $q_i$. Its acceptance probability is positive, so it terminates almost surely, although it can be inefficient.

An efficient exact alternative is adaptive rejection sampling. The log density hasThus it is a log-concave probability density. Choose one tangent point with positive derivative in the left tail and one with negative derivative in the right tail; these exist because the Gaussian term dominates there. Their tangent upper hull gives an integrable piecewise exponential rejection envelope. Sample that envelope, accept using the target-to-envelope ratio, and add evaluated points to improve it. Concavity guarantees the upper bound, so accepted samples remain exact. Sample each coordinate independently by either method; a finite number of Metropolis steps would not provide the requested exact draws.