We introduce PHI, a fully Bayesian Markov chain Monte Carlo algorithm designed for the structural decomposition of galaxy images. PHI uses a triple layer approach to effectively and efficiently explore the complex parameter space. Combining this with the use of priors to prevent non-physical models, PHI offers a number of significant advantages for estimating surface brightness profile parameters over traditional optimization algorithms. We apply PHI to a sample of synthetic galaxies with Sloan Digital Sky Survey (SDSS)-like image properties to investigate the effect of galaxy properties on our ability to recover unbiased and well-constrained structural parameters. In two-component bulge+disc galaxies, we find that the bulge structural parameters are recovered less well than those of the disc, particularly when the bulge contributes a lower fraction to the luminosity, or is barely resolved with respect to the pixel scale or point spread function (PSF). There are few systematic biases, apart from for bulge+disc galaxies with large bulge Sérsic parameter, n. On application to SDSS images, we find good agreement with other codes, when run on the same images with the same masks, weights, and PSF. Again, we find that bulge parameters are the most difficult to constrain robustly. Finally, we explore the use of a Bayesian information criterion method for deciding whether a galaxy has one or two components.