In this thesis I develop Bayesian approach to model galaxy surface brightness and apply it to a bulge-disc decomposition analysis of galaxies in near-infrared band, from Two Micron All Sky Survey (2MASS). The thesis has three main parts. First part is a technical development of Bayesian galaxy image decomposition package Galphat based on Markov chain Monte Carlo algorithm. I implement a fast and accurate galaxy model image generation algorithm to reduce computation time and make Bayesian approach feasible for real science analysis using large ensemble of galaxies. I perform a benchmark test of Galphat and demonstrate significant improvement in parameter estimation with a correct statistical confidence. Second part is a performance test for full Bayesian application to galaxy bulgedisc decomposition analysis including not only the parameter estimation but also the model comparison to classify different galaxy population. The test demonstrates that Galphat has enough statistical power to make a reliable model inference using galaxy photometric survey data. Bayesian prior update is also tested for parameter estimation and Bayes factor model comparison and it shows that informative prior significantly improves the model inference in every aspects. Last part is a Bayesian bulge-disc decomposition analysis using 2MASS Ks-band selected samples. I characterise the luminosity distributions in spheroids, bulges and discs separately in the local Universe and study the galaxy morphology correlation, by full utilising the ensemble parameter posterior of the entire galaxy samples. It shows that to avoid a biased inference, the parameter covariance and model degeneracy has to be carefully characterised by the full probability distribution.