Past exam of the mathematics course of the University of Cambridge / 2019 / iii / Paper 219 (Astrostatistics)

LetSubtracting the measured distance modulus from the measured apparent magnitude givesFor the Hubble flow, defineThe Hubble law, with the Hubble constant parametrized by $\theta =5{\log }_{10}h$, gives $r_i\sim N(M_0-\theta ,H_i)$. After integrating over each intrinsic absolute magnitude and each unobserved true distance modulus, independence therefore gives the likelihood function

Put $w_k=C_k^{-1}$, $u_i=H_i^{-1}$, $A={\sum }_k{w}_k$, $B={\sum }_i{u}_i$, and form the weighted meansThe two equations obtained from the score function areConsequently the maximum-likelihood estimators areThe Hessian matrix of the log likelihood isIts first leading principal minor is negative and its determinant is $AB>0$, so it is a negative-definite matrix. Thus the stationary point is the unique global maximum.

Under homoskedasticity, writeThen ${\hat{M}}_0=\bar{q}$ and $\hat{\theta }=\bar{q}-\bar{r}$. Both are unbiased estimators, and their covariance matrix isIndeed the Fisher information isand its inverse is exactly the displayed covariance matrix. The estimators therefore attain the multivariate Cramer-Rao bound and are efficient estimators.

By the invariance property of maximum likelihood estimation,The sampling distribution is $\hat{\theta }\sim N(\theta ,{\sigma }_{\hat{\theta }}^2)$, where$\hat{h}$ has a log-normal distribution withWriting $v={\sigma }_{\hat{\theta }}^2/{\alpha }^2$, its exact fractional bias and variance areHence, to leading order,

The measurement model has the Gaussian likelihoodMarginalizing the latent angular momentum gives the normalized posterior distributionIts posterior mean is

Treat the simulated pairs $(x_i,m_i)$ as samples from the prior $P(x,m)$. The numerator and denominator of the posterior mean are then ordinary Monte Carlo estimators, sowhere the normalized importance sampling weights areThe common Gaussian normalizing constant cancels.

For arbitrary nonnegative raw weights $W_i$, let $\bar{W}=K^{-1}{\sum }_i{W}_i$. The empirical squared coefficient of variation, using variance divisor $K$, isSubstitution into the stated definition gives the usual effective sample size of importance samplingFor normalized weights $w_i=W_i/{\sum }_j{W}_j$, this reduces to $\boxed{\widehat{\mathrm{ESS}}=1/\sum _i{w}_i^2}$.

Conditional independence givesFor each satellite, Bayes theorem gives $P(d_s\mid m)\propto P(m\mid d_s)/P(m)$. ThereforeIf ${\hat{p}}_0(m)$ is a kernel density estimator for the simulated marginal masses and ${\hat{p}}_s(m)$ estimates the posterior based on satellite $s$, thenThe one-dimensional integrals can be evaluated by numerical integration on a common mass grid.

Let $p(m)=P(m\mid \mathbf{d})$ and draw independently from an importance density $q$. The unbiased estimatorof $I=\int mp(m)dm$ has one-sample second momentBy the Cauchy-Schwarz inequality,Equality holds precisely when $q(m)\propto |m|p(m)$, giving the optimal importance density for a single integralThis is circular in practice: constructing and normalizing $q^{\ast }$ requires detailed knowledge of the posterior and the expectation of $|m|$. Here log masses are positive, so the unknown normalizer is the posterior mean being estimated. It is also optimal only for this one integral, not for general posterior summaries.

Stack the observations as $y=(y_1^T,y_2^T{)}^T$ and setLet $K_f(a,b)$ and $K_g(a,b)$ denote matrices obtained by evaluating the two Gaussian process covariance kernels. Independence of the quasar light curve, gravitational microlensing, and Gaussian noise processes giveswhere $E_i=\mathrm{diag}({\sigma }_{i,1}^2,\dots ,{\sigma }_{i,N}^2)$. Thus $y\mid \theta$ is a multivariate normal distribution and its Gaussian-process marginal likelihood isThe off-diagonal blocks are essential: both images contain the same delayed Ornstein-Uhlenbeck process.

At the fitted parameters, let $C=C_{\hat{\theta }}$ and $\mu ={\mu }_{\hat{\theta }}$. For prediction times $t_{\ast }$ defineThe microlensing processes and measurement errors contribute no cross-covariance with the latent quasar light curve. The Gaussian process regression posterior is thereforeThe requested pointwise posterior variances are the diagonal entries of the latter matrix.

Use broad proper uniform priors for $\Delta t$, $\Delta m$, and $c$ over physically plausible ranges, and broad log-uniform priors for the positive scales $A_f$ and ${\tau }_f$. ThenA Random-walk Metropolis algorithm can update $(\Delta t,\Delta m,c,\log A_f,\log {\tau }_f)$ with a multivariate Gaussian proposal distribution. Initialize several dispersed chains near plausible cross-correlation delays and near the marginal-likelihood optimum; reject proposals outside the prior bounds; discard warm-up while adapting only the proposal scale and covariance; then freeze the kernel and retain a long run. Evaluate trace plots, autocorrelations, acceptance rates, between-chain agreement, and the effective sample size of a Markov chain. Posterior predictive quasar light curves provide a model check.

Write the target posterior density as $\pi (\theta )$ and the proposal distribution density as $q({\theta }^{\prime }\mid \theta )$. The Metropolis–Hastings algorithm accepts a proposed move withFor distinct states,which is symmetric in $\theta$ and ${\theta }^{\prime }$. The rejection probability supplies the diagonal part, so the entire transition kernel satisfies detailed balance. Integrating the detailed-balance identity over the starting state proves $\int \pi (\theta )P(\theta ,d{\theta }^{\prime })=\pi ({\theta }^{\prime })d{\theta }^{\prime }$. Hence the posterior is a stationary distribution; an irreducible Markov chain that is also an aperiodic Markov chain converges uniquely to it.

For one object, the probabilistic graphical model factorization isEach factor is the normal distribution density specified by the model. This factorization displays the conditional independences of the latent variables ${\xi }_i,{\eta }_i$ and the noisy observations $x_i,y_i$.

With the stated flat priors, the full joint density, up to a constant, isThe priors on $\alpha ,\beta ,$ and $\mu$ contribute constants on $\mathbb{R}$, while those on the two variances contribute constants on $(0,\infty )$. These are improper priors, so posterior propriety must be checked; the full-rank, sufficiently large-data case used below is proper.

For each $i$, place ${\xi }_i$ and ${\eta }_i$ inside a plate replicated $N$ times. The directed edges are represented byThe shaded observed nodes are $x_i,y_i$; the unshaded nodes ${\xi }_i,{\eta }_i$ are latent; and $\alpha ,\beta ,{\sigma }^2,\mu ,{\tau }^2$ lie outside the plate. This is the probabilistic graphical model encoded by the joint factorization.

Every move can be drawn from a full conditional distribution, producing a Gibbs sampler with acceptance probability one. Write $s_{x,i}^2={\sigma }_{x,i}^2$ and $s_{y,i}^2={\sigma }_{y,i}^2$. First update independentlyand thenLet $X$ have rows $(1,{\xi }_i)$ and $\eta =({\eta }_i)$. Update the linear regression coefficients jointly byand updateFinally, the flat positive variance priors give the following full conditionals, each an inverse-gamma distribution:A systematic sweep in the displayed order, using the newly sampled values immediately, defines the chain. The shapes are positive for $N>2$; full column rank of $X$ is also required.