We present a Bayesian method for inferring axisymmetric plasma equilibria from the magnetic field and plasma pressure measurements. The method calculates all possible solutions for plasma current and pressure distributions consistent with the measurements and magnetohydrodynamic (MHD) force balance. Toroidal plasma current and magnetic field coils are modelled as a set of axisymmetric current-carrying solid beams. The other parameters such as plasma pressure and poloidal current flux are given as a function of poloidal magnetic flux, which is determined given a 2D current distribution. Plasma pressure and poloidal current flux profiles are modelled as Gaussian processes whose smoothness is optimally chosen based on the principle of Occam's razor. To find equilibrium solutions, we introduce an MHD force balance constraint at every plasma current beam as a part of the prior knowledge. Given all these physical quantities, predictions calculated by the predictive (forward) models for diagnostics are compared to the observations. The high dimensional complex posterior probability distribution is explored by a new algorithm based on the Gibbs sampling scheme.