For each quasar image, we run the MCMC with 200 chains, and a minimum of 10,000 iterations with the first 5,000 discarded as a burn-in period (systems with more surrounding galaxies required up to twice as many iterations to obtain convergence). To ensure that the model is well-fit to the data, we examine the resulting posterior distributions, altering the allowed parameter range and iteration count until each parameter has converged and the residual flux in the model subtracted image is consistent with random noise.

For the six quasars and their companion galaxies, we create posterior-weighted model images before convolution with the PSF, and after the model has been convolved with the PSF and subtracted from the original image. These weighted images are the (per-pixel) mean of all sample images, with more probable locations in parameter space being more densely populated with samples. The resulting images for the six quasars in the J- and H-bands are shown in Figures \[fig:quasar1\_resids\]–\[fig:quasar3\_resids\]. The residual images show a central core of flux, which contains some residual flux from the quasar, and may also contain underlying host emission. The regions in the residual images where companion galaxy models have been subtracted purely consists of noise, 