Stack the observations as and set
Let and denote matrices obtained by evaluating the two Gaussian process covariance kernels. Independence of the quasar light curve, gravitational microlensing, and Gaussian noise processes gives
where . Thus is a multivariate normal distribution and its Gaussian-process marginal likelihood is
The off-diagonal blocks are essential: both images contain the same delayed Ornstein-Uhlenbeck process.
Use broad proper uniform priors for , , and over physically plausible ranges, and broad log-uniform priors for the positive scales and . Then
A Random-walk Metropolis algorithm can update 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.