Solution (source code)

= Solution

Write $v=\sigma_{\rm int}^2$, $r_s=\sigma_{m,s}^2$ and $d_s=\mu(z_s;H_0,w,\Omega_M)$. The latent <absolute magnitude> and measurement error are independent normal variables, so their <convolution> is normal. After integrating out the individual intrinsic magnitudes, each <apparent magnitude> has distribution
$$
\widehat m_s\mid z_s,M_0,v,H_0,w,\Omega_M\sim N(M_0+d_s,V_s),\qquad V_s=v+r_s.
$$
Consequently the <likelihood function> for the independent sample is
$$
\boxed{L=\prod_{s=1}^{N}(2\pi V_s)^{-1/2}\exp\left[-\frac{(\widehat m_s-M_0-d_s)^2}{2V_s}\right]}.
$$
This is the marginal observational likelihood of a <hierarchical Bayesian model>, not the likelihood conditional on each unknown intrinsic magnitude. The individual errors remain <heteroscedastic> through their known $r_s$.