  We describe an efficient and exact method that enables global Bayesian
analysis of cosmic microwave background (CMB) data. The method reveals the
joint posterior density (or likelihood for flat priors) of the power spectrum
$C_\ell$ and the CMB signal. Foregrounds and instrumental parameters can be
simultaneously inferred from the data. The method allows the specification of a
wide range of foreground priors. We explicitly show how to propagate the
non-Gaussian dependency structure of the $C_\ell$ posterior through to the
posterior density of the parameters. If desired, the analysis can be coupled to
theoretical (cosmological) priors and can yield the posterior density of
cosmological parameter estimates directly from the time-ordered data. The
method does not hinge on special assumptions about the survey geometry or noise
properties, etc. It is based on a Monte Carlo approach and hence parallelizes
trivially. No trace or determinant evaluations are necessary. The feasibility
of this approach rests on the ability to solve the systems of linear equations
which arise. These are of the same size and computational complexity as the
map-making equations. We describe a pre-conditioned conjugate gradient
technique that solves this problem and demonstrate in a numerical example that
the computational time required for each Monte Carlo sample scales as
$n_p^{3/2}$ with the number of pixels $n_p$. We test our method using the
COBE-DMR data and explore the non-Gaussian joint posterior density of the
COBE-DMR $C_\ell$ in several projections.
