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

An apparent magnitude is the sum of absolute magnitude and distance modulus. ThusThe independent unbiased distance-modulus estimate is $N({\mu }_{\mathrm{LMC}},{\sigma }_{\mathrm{LMC}}^2)$, so closure of independent normal distributions under linear combinations givesHence $\mathbb{E}{\hat{S}}_1=M_0^{\ast }$ and $\mathrm{Var}({\hat{S}}_1)=V_1$.

The common distance modulus of calibrator galaxy $g$ cancels between its supernova and mean Cepheid apparent magnitudes. Their independent intrinsic scatters therefore giveThus $\mathbb{E}{\hat{S}}_2^g=\Delta M$ and $\mathrm{Var}({\hat{S}}_2^g)=V_2$.

The Hubble law gives $d_i=c{z}_i/(100h)$ Mpc, so its distance modulus isConsequently

All three sets of statistics are independent. With $s_1={\hat{S}}_1$, $s_{2g}={\hat{S}}_2^g$, and $s_{3i}={\hat{S}}_3^i$, their likelihood function, up to a parameter-independent factor, is

Differentiating the log-likelihood gives the unique stationary pointIts Hessian is diagonal with entries $-1/V_1$, $-G/V_2$, and $-N_{\mathrm{SN}}/{\sigma }_{\mathrm{SN}}^2$, so it is the unique maximum. The estimators are unbiased and have varianceswhich equal their Cramér-Rao lower bounds because the normal location statistics are efficient.

Since $\theta =M_0^{\ast }+\Delta M-\mathcal{M}$, invariance of the maximum-likelihood estimator givesIt is unbiased and normal with variance

By invariance of maximum likelihood,Because $\hat{\theta }-\theta \sim N(0,{\sigma }_{\theta }^2)$, the ratio $\hat{h}/h$ is log-normal. Its exact variance isAs $N_{\mathrm{SN}},G,N_{\mathrm{LMC}}\to \infty$, all sampling terms vanish but $V_1$ retains ${\sigma }_{\mathrm{LMC}}^2$. The dominant remaining uncertainty is therefore the parallax calibration of the LMC distance modulus.

Writing $\phi (u;m,v)$ for the $N(m,v)$ density, the conditional-independence structure of the latent-variable model gives

First integrate out ${\eta }_i$: conditional on ${\xi }_i$, $y_i\sim N(\alpha +\beta {\xi }_i,s^2)$ with $s^2={\sigma }^2+{\sigma }_y^2$. Since $(x_i,y_i)$ is an affine transformation of independent normal variables, it is bivariate normal with meanand covarianceIts determinant simplifies toFor $r_i=(x_i-\mu ,y_i-\alpha -\beta \mu {)}^T$, the observed-data likelihood is therefore

Take flat priors on $\alpha ,\beta ,\mu$ and scale priors $\pi ({\sigma }^2,{\tau }^2)\propto ({\sigma }^2{\tau }^2{)}^{-1}$ on the positive variances. The posterior distribution is thenA random-walk Metropolis–Hastings algorithm can update $(\alpha ,\beta ,\mu ,\log {\sigma }^2,\log {\tau }^2)$ with a symmetric proposal and accept a proposed state $z^{\prime }$ from $z$ with probability $1\wedge \pi (z^{\prime }\mid \mathcal{D})/\pi (z\mid \mathcal{D})$, including the Jacobian if the target is represented in transformed coordinates. Its transition kernel satisfieswhich is detailed balance; hence the posterior is stationary.

The marginal predictor density isGaussian conditional expectation givesIntegrating this conditional distribution through the response model yields

When $\tau \gg {\sigma }_x$, $\kappa \to 1$, $v\to {\sigma }_x^2$, and $A\sim {\tau }^2$. HencewhereThe approximate maximum-likelihood estimators are

Put $\rho =e^{-\Delta t/\tau }=e^{-x}$ and retain $R=k(0,0)$, which equals one for the stated exponential covariance function. The covariance matrix of $(y_1,y_2,y_3)$ isConditioning a multivariate normal distribution and simplifying givesThe absence of $y_1$ from both conditional moments shows that $y_3$ and $y_1$ are conditionally independent given $y_2$. Equivalently, the exponential-kernel Gaussian process is the stationary Ornstein-Uhlenbeck process, which is Markov.

As $\Delta t/\tau \to \infty$, $\rho \to 0$. Therefore the posterior predictive mean tends to $\mu$ and its variance tends to $R$: a sufficiently distant observation has reverted to the stationary prior distribution.

The Markov factorization isDifferentiating its log-likelihood with respect to $\mu$ givesIts coefficients sum to one, so it is unbiased. Direct covariance calculation, equivalently inversion of its Fisher information, givesThus the variance tends to $R$ as $\Delta t/\tau \to 0$, because the observations become perfectly correlated, and to $R/3$ as $\Delta t/\tau \to \infty$, because they become three independent draws.

One simple algorithm is importance sampling from the prior. Draw $({\theta }_j,{\mathbf{\alpha }}_j)$ independently from $\pi (\theta )\pi (\mathbf{\alpha })$ and assign weight $w_j=L_A({\theta }_j,{\mathbf{\alpha }}_j)$. Thenestimates the Bayesian model evidence, and the normalized weights $w_j/{\sum }_k{w}_k$ represent the posterior. A weighted histogram or kernel density estimation of the ${\theta }_j$ estimates the marginal ${\mathcal{P}}_A(\theta )$; discarding ${\mathbf{\alpha }}_j$ performs the marginalization.

Independence of the experiments giveswhereThe marginal posterior is obtained by integrating the displayed joint posterior over both nuisance-parameter vectors.

Define the prior-predictive nuisance integralsThenFor each individual analysis, Bayes theorem also gives

Estimate each one-dimensional marginal posterior from its individual experiment's weighted samples, for example by weighted kernel density estimation. Combining that density estimate with the individual evidence estimate givesTheir product with $\pi (\theta )$ can be normalized on the one-dimensional $\theta$ space, avoiding all joint nuisance-parameter sampling.

Compute the joint evidence by one-dimensional numerical integration,Equivalently, the individual outputs giveQuadrature or one-dimensional importance sampling evaluates this integral without entering the $(\theta ,\mathbf{\alpha },\mathbf{\beta })$ space.