  We present a general method for Bayesian inference of the underlying
covariance structure of random fields on a sphere. We employ the Bipolar
Spherical Harmonic (BipoSH) representation of general covariance structure on
the sphere. We illustrate the efficacy of the method as a principled approach
to assess violation of statistical isotropy (SI) in the sky maps of Cosmic
Microwave Background (CMB) fluctuations. SI violation in observed CMB maps
arise due to known physical effects such as Doppler boost and weak lensing; yet
unknown theoretical possibilities like cosmic topology and subtle violations of
the cosmological principle, as well as, expected observational artefacts of
scanning the sky with a non-circular beam, masking, foreground residuals,
anisotropic noise, etc. We explicitly demonstrate the recovery of the input SI
violation signals with their full statistics in simulated CMB maps. Our
formalism easily adapts to exploring parametric physical models with non-SI
covariance, as we illustrate for the inference of the parameters of a Doppler
boosted sky map. Our approach promises to provide a robust quantitative
evaluation of the evidence for SI violation related anomalies in the CMB sky by
estimating the BipoSH spectra along with their complete posterior.
