Bayesian fitting of multi-Gaussian expansion models to galaxy images

Fitting parameterized models to images of galaxies has become the standard for measuring galaxy morphology. This forward modelling technique allows one to account for the PSF to effectively study semi-resolved galaxies. However, using a specific parameterization for a galaxy's surface brightness profile can bias measurements if it is not an accurate representation. Furthermore, it can be difficult to assess systematic errors in parameterized profiles. To overcome these issues we employ the Multi-Gaussian expansion (MGE) method of representing a galaxy's profile together with a Bayesian framework for fitting images. MGE flexibly represents a galaxy's profile using a series of Gaussians. We introduce a novel Bayesian inference approach which uses pre-rendered Gaussian components, which greatly speeds up computation time and makes it feasible to run the fitting code on large samples of galaxies. We demonstrate our method with a series of validation tests. By injecting galaxies, with properties similar to those observed at $z\sim1.5$, into deep HST observations we show that it can accurately recover total fluxes and effective radii of realistic galaxies. Additionally we use degraded images of local galaxies to show that our method can recover realistic galaxy surface brightness and color profiles. Our implementation is available in an open source python package $\texttt{imcascade}$, which contains all methods needed for the preparation of images, fitting and analysis of results.