Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 216 (Bayesian Modelling and Computation)

Let $P$ be the Markov kernel and $\mu$ its invariant probability measure. Using the conventiongeometric ergodicity means that there are $r\in (0,1)$ and a finite measurable state-dependent prefactor $M(x)$ such thatHere total variation distance is used, and the rate $r$ does not depend on $x$. The prefactor need not be uniformly bounded; that stronger property is uniform geometric ergodicity. For a deterministic starting state, this estimate must hold at that state. If the definition is made only $\mu$-almost everywhere, the subsequent assertions concern those starting states where the prefactor is finite. The reversible Markov chain assumption means $\mu (dx)P(x,dy)=\mu (dy)P(y,dx)$.

Put $B=\Vert f{\Vert }_{\infty }$, $v={\mathrm{Var}}_{\mu }(f)$, and $N_n=n-m_n$. We first obtain a uniform-in-$f$ stationary bound, then compare the burned-in trajectory with a stationary one.

For a reversible Markov chain, geometric ergodicity gives an absolute L2 spectral gap of a reversible Markov chain. Thus, for some $r<1$, the Markov kernel viewed as an operator satisfiesOne can use the rate in the total-variation estimate. To see why, reversibility makes $P$ a self-adjoint operator. For bounded mean-zero $h$ supported on a set where $M$ is bounded,The spectral theorem for normal operators writes the left side as $\int t^{2k}d{\nu }_h(t)$. Any positive spectral mass where $|t|>r$ would contradict the inequality as $k\to \infty$. These test functions are dense in $L_0^2(\mu )$: truncate a function and restrict it to $\{M\le j\}$, then correct its mean using a bounded function supported on a fixed positive-measure set where $M$ is bounded. Hence the spectral restriction extends to the entire mean-zero space.

Let $(W_i)$ be the stationary chain. For $h=f-\mu f$, the stationary covariance bound for reversible Markov chains givesConsequently

The initial law of the burned-in trajectory is $P^{m_n}(x,\cdot )$. Its total variation distance from $\mu$ is at most ${\epsilon }_n=M(x)r^{m_n}$. Attaching the same conditional future evolution cannot increase that distance, so the entire length-$N_n$ path laws differ by at most ${\epsilon }_n$. For either path the sum $S$ has $|S|\le B{N}_n$. The dual bound for total variation distance gives a difference of at most $2B^2{N}_n^2{\epsilon }_n$ in the second moments, and at most $4B^2{N}_n^2{\epsilon }_n$ in the squares of the means. This burn-in total variation comparison yieldsSince $m_n=\lfloor n^{1/3}\rfloor$, $N_n{r}^{m_n}\to 0$. Dividing by $N_n$ givesThe constant depends on the chain, not on $f$.

Let $A_n={\sum }_{i=1}^{m_n}f(X_i)$ and $T_n={\sum }_{i=m_n+1}^{n}f(X_i)$. Since $|A_n|\le B{m}_n$, $\mathrm{Var}(A_n)\le B^2{m}_n^2$. The Cauchy-Schwarz inequality applied to the centered sums givesDivide by $\sqrt{n}$. The first term is at most $B{m}_n/\sqrt{n}\to 0$, while part (b) and $(n-m_n)/n\to 1$ give a limsup of at most $\sqrt{\gamma {\mathrm{Var}}_{\mu }(f)}$ for the second. Squaring provesThis argument also covers ${\mathrm{Var}}_{\mu }(f)=0$; it does not require independence between the discarded prefix and the retained path. The same chain-dependent $\gamma$ as in part (b) suffices.

For the planar uniform distribution with a Lebesgue measure density, assume that the compact convex set $\mathcal{C}$ has nonempty interior, equivalently positive area. This assumption is missing from the wording: a line segment or singleton would not have such a planar uniform density.

The hit-and-run sampler at a current point $x$ chooses $u=(\cos \Theta ,\sin \Theta )$ with $\Theta$ uniform on $[0,2\pi )$, finds the intervaland samples $T$ uniformly on that interval, setting $x^{\prime }=x+Tu$. Convexity makes the intersection an interval and compactness makes it bounded.

If $x$ is in the interior, every direction has positive chord length, and the next point is in the interior almost surely. For a definition valid at every boundary point as well, redraw the direction whenever the chord has zero length. We use this convention in the following parts. Positive-length directions have positive angular probability at every $x$, since an interior ball is visible from $x$, so the redraw terminates almost surely.

The Markov kernel has the uniform distribution on $\mathcal{C}$ as its invariant law, as follows from the density symmetry proved next. Repeating the step therefore produces approximate uniform samples after convergence. The boundary convention matters: if singleton chords were instead treated as holding moves, a square's corner would have a positive atom at its starting point, contradicting part (b)'s claim of an everywhere absolutely continuous kernel.

Write $\ell (x,u)=t_{+}(x,u)-t_{-}(x,u)$ and letThis is the accepted fraction of directions when zero-length chords are redrawn; it equals one in the interior. The Markov kernel integrates a test function $h$ asThe signed polar representation has Jacobian determinant $|t|$. Each point $y\ne x$ is represented twice, once with positive $t$ and once with negative $t$ and the opposite direction. Since the chord is the same in either representation, the hit-and-run kernel density isfor almost every $y\in \mathcal{C}$, and it is zero outside the body. Thus the transition is absolutely continuous with respect to planar Lebesgue measure. Redrawing removes the only potential zero-chord atoms.

Let $D=\mathrm{diam}(\mathcal{C})>0$. For distinct $x,y\in \mathcal{C}$, the chord contains their segment, so it has positive length, andHence $p(x,y)\ge 1/(\pi D^2)$. Values on the diagonal, and on other null sets where a density is only specified almost everywhere, can be assigned this same lower bound. This gives a version satisfyingFor interior $x,y$, $a(x)=a(y)=1$ and the chord length is unchanged on interchanging its endpoints, so $p(x,y)=p(y,x)$. A convex body's boundary has area zero; for example contractions about an interior point fill the interior and their areas tend to the body's area. Thus this symmetry holds almost everywhere for two uniform points, proving detailed balance and stationarity of the uniform distribution. Without positive area or the boundary convention in part (a), the literal assertion need not hold.

The geometric drift condition requires a measurable $V:\mathcal{X}\to [1,\infty )$, a small set $C$, constants $\lambda \in (0,1)$ and $b<\infty$, andSmallness means that for some $m\ge 1$, $\epsilon >0$, and probability measure $\nu$, $P^m(x,\cdot )\ge \epsilon \nu (\cdot )$ for every $x\in C$. With irreducibility and aperiodicity, this drift-minorisation condition implies geometric ergodicity. The function $V$ controls the chain away from the small set; inside that set the additive constant allows unrestricted bounded drift.

Let $A=\mathrm{Leb}(\mathcal{C})>0$ and let $\mu$ be the uniform distribution on the body. The hit-and-run kernel density lower bound impliesFor the last range, $\mathcal{C}$ lies in a radius-$D$ disc about any of its points, so $A\le \pi D^2$ and $\epsilon \le 1/2$. This global Doeblin condition proves irreducibility with respect to $\mu$ and aperiodicity. The entire state space is a small set; $V=1$, $\lambda =1/2$, and $b=1/2$ satisfy the geometric drift condition.

In fact the conclusion is stronger. Write $P=\epsilon \Pi +(1-\epsilon )R$, where $\Pi (x,\cdot )=\mu$ and $R$ is the residual Markov kernel. Since $\mu P=\mu$, also $\mu R=\mu$. Each transition has a probability $\epsilon$ of resetting to stationarity, and after such a reset every subsequent distribution remains stationary. Equivalently,Thus the algorithm has uniform geometric ergodicity, using the full-dimensional and boundary conventions stated in part (a).

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.

In factor analysis, the independent multivariate normal distributions of $Z_i$ and ${\xi }_i$ imply that their linear combination is also multivariate normal. Its expected value is zero and its covariance matrix isThereforeThe noise variance is assumed positive. The factor contribution is positive semidefinite and has rank at most $k$, while the noise makes the marginal covariance nonsingular.

For one factor, write $\lambda \in {\mathbb{R}}^p$ for the loading vector and $v={\sigma }^2>0$. Let ${\lambda }_t$ be its current value and $s=v+\Vert {\lambda }_t{\Vert }^2$. The E-step of the expectation-maximization algorithm uses the conditional multivariate normal distribution to obtainThus $\mathbb{E}[Z_i^2\mid Y_i]=u+m_i^2$. The terms of the expected complete-data log-likelihood that depend on a candidate $\lambda$ areThis is a strictly concave quadratic in $\lambda$, since ${\sum }_i(u+m_i^2)>0$. Its unique maximizer isFor an explicit expression involving only the current loading and data, put $S=n^{-1}{\sum }_i{Y}_i{Y}_i^T$. Substitution gives the EM update for a single-factor Gaussian modelNo empirical centering is applied because the model fixes the mean to zero. If initialized at zero, this update remains at zero; monotonicity alone does not guarantee that such a stationary point is a global maximum. This part is present in the PDF but missing from the TeX transcription.

Fix the observed $Y$, let $L(\theta )=p_{\theta }(Y)$, and let $q_0(z)=p_{{\theta }_0}(z\mid Y)$. The E-step formsFor finite expected log-likelihoods and $L({\theta }_0)>0$, the Jensen inequality givesThe first inequality allows additional mass outside the old complete-data support; it is equality when the relevant supports coincide. Ratios are taken on the support of $q_0$. If a candidate assigns zero density on positive $q_0$-mass, its $Q$ is minus infinity and it cannot be an improving finite M-step.

An M-step of the expectation-maximization algorithm chooses ${\theta }_1$ with $Q({\theta }_1\mid {\theta }_0)\ge Q({\theta }_0\mid {\theta }_0)$, soThis EM likelihood monotonicity also holds for a generalized M-step that merely increases $Q$, rather than maximizing it exactly. It is a statement about the observed-data likelihood function.

Let ${\varphi }_d(z)=(2\pi {)}^{-d/2}\exp (-\Vert z{\Vert }^2/2)$ be the standard multivariate normal density. Independence gives the proposal distribution densityFor an orthogonal matrix $A$,The isotropic Gaussian density is unchanged by this transformation. Since $A^T$ has the same distribution as $A$,Thus this is an orthogonal-mixture Metropolis proposal with a symmetric density; inversion invariance, rather than any assumption of a uniform distribution on orthogonal matrices, is what is needed.

With $\alpha (x,y)=\min \{1,\pi (y)/\pi (x)\}$, the accepted off-diagonal transitions satisfyThe rejection part is supported on the diagonal and satisfies detailed balance automatically. Values of the acceptance rule at starting points where $\pi (x)=0$ can be specified separately; they do not affect stationarity under $\pi$. Therefore the Metropolis–Hastings algorithm kernel is reversible with invariant density $\pi$, proving

Let $h(\beta )=\Phi (x_{\mathrm{test}}^T\beta )$. To distinguish the number of posterior draws from the training-sample size, write it as $N$; in the question's sample notation $N=n$. The Monte Carlo estimator of the requested posterior mean isFor independent exact posterior draws it is unbiased, andThe features are conditioned on as fixed data. The estimated probability concerns a future genuine banknote under the posterior predictive probability, not an indicator sampled from that predictive distribution.

Put $t_i=x_i^T\beta$ and $s_i=2Y_i-1\in \{-1,1\}$. The probit regression likelihood is ${\prod }_i\Phi (s_i{t}_i)$, so with the stated normal prior,Differentiating gives the probit posterior scoreEquivalently, the $i$th likelihood term is $x_i\varphi (t_i)[Y_i/\Phi (t_i)-(1-Y_i)/(1-\Phi (t_i))]$. The prior's gradient is essential.

This posterior density is smooth, positive everywhere, and bounded above by a constant times the Gaussian prior, since the likelihood is at most one and has positive normalizing constant. Its derivatives are integrable: differentiating a likelihood factor gives a bounded normal density, and differentiating the prior gives a linear factor times a Gaussian. Integration by parts therefore yieldsThus each component has a known zero expected value and can be used as a posterior score control variate. This identity differentiates the posterior density in the random parameter; it is distinct from the usual mean-zero score identity for a likelihood under repeated sampling.

The requisite second moments are finite. For $t\ge 0$, $\varphi (t)/\Phi (t)$ is bounded. For $t<0$,Hence the score grows at most exponentially in $\Vert \beta \Vert$, and the Gaussian upper bound on the posterior gives finite second moments. For a fixed vector $a$, averaging $h(\beta )-a^{T}g(\beta )$ is therefore an unbiased Monte Carlo estimator of the same mean; suitable coefficients can reduce its variance.

All moments below are under the proper posterior distribution of part (b). Write the known covariance matrix in block formThe score covariance $S$ is positive definite in this model. Indeed, a further integration by parts gives $\mathbb{E}[g(\beta ){\beta }^T]=-I_p$. If $a^{T}Sa=0$, then $a^{T}g=0$ almost surely, because $\mathbb{E}[g]=0$. Multiplying by $a^T\beta$ would give $0=-\Vert a{\Vert }^2$, so $a=0$.

For a fixed coefficient vector $a$, the control variate summand has varianceThus the unique minimum is attained at $a_{\ast }=S^{-1}c$, givingThe coefficient is fixed because the covariance is assumed known, and the posterior score control variate has mean zero, so this estimator remains unbiased.

We can identify exactly when the improvement is strict. Integration by parts with the bounded smooth function $h$ givesThe scalar expectation is strictly positive. Therefore $c\ne 0$ whenever $x_{\mathrm{test}}\ne 0$, and positive definiteness gives a strictly smaller variance in that case. If $x_{\mathrm{test}}=0$, $h\equiv 1/2$ and both estimators already have variance zero. The printed request for a smaller variance therefore needs this nondegeneracy qualification; a non-increasing variance always holds.