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 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$.

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.

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.

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.

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.

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.

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.

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}$.

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.

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

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,

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.

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.

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

The MA(1) model has much smaller AIC, $2871.79$ versus $3217.1$, so select MA(1). A nominal Wald 95% interval for its non-intercept parameter isThe estimate is close to the noninvertible boundary $\theta =1$, where the regular asymptotic normal approximation becomes poor and likelihood curvature can understate the true one-sided uncertainty. Therefore statement (iii) is the most plausible: the nominal interval is too narrow to attain its stated coverage reliably.

An MA($q$) process is invertible when its innovations admit a causal absolutely summable linear representation in present and past observations. With the backshift operator $B$, MA(1) satisfiesIf $|\theta |<1$, the geometric series converges absolutely:Thus the process is invertible.

An MA($q$) process iswhere $({\epsilon }_t)$ is white noise of variance ${\sigma }^2$. For MA(1), $X_t={\epsilon }_t+\theta {\epsilon }_{t-1}$ andThe transformation $(\theta ,{\sigma }^2)\mapsto (1/\theta ,{\theta }^2{\sigma }^2)$ leaves this autocovariance unchanged. A Gaussian process is determined by its mean and covariance, so the parameters are not identifiable unless one selects, for example, the invertible representative $|\theta |<1$.