Muon tomography is an emerging imaging technique using muons contained in natural cosmic rays to determine subsurface density in a diverse range of objects. It is often difficult to get that information from other imaging techniques and muon tomography can be carried out in a cost-effective and environmentally sound manner. The process of inference from this muon tomography data can be framed as an inverse problem, usually with non-unique solutions. In this paper, we carry out statistical inversion in a Bayesian regression framework for the subsurface density contrast based on data from muon detectors, first in the case where there is a single uniform anomaly, and extend it to multiple anomalies. Based on a Gaussian prior model we have a closed form solution for the posterior distribution that can be efficiently implemented. A posterior uncertainty representation has been developed directly from the Bayesian inverse solution which gives an uncertainty quantification for the computed solution. We also compare our Bayesian inversion approach with industry standard regularization schemes. The method is illustrated in a series of synthetic examples and a real-world application.