  To accomplish correct Bayesian inference from weak lensing shear data
requires a complete statistical description of the data. The natural framework
to do this is a Bayesian Hierarchical Model, which divides the chain of
reasoning into component steps. Starting with a catalogue of shear estimates in
tomographic bins, we build a model that allows us to sample simultaneously from
the the underlying tomographic shear fields and the relevant power spectra
(E-mode, B-mode, and E-B, for auto- and cross-power spectra). The procedure
deals easily with masked data and intrinsic alignments. Using Gibbs sampling
and messenger fields, we show with simulated data that the large (over
67000-)dimensional parameter space can be efficiently sampled and the full
joint posterior probability density function for the parameters can feasibly be
obtained. The method correctly recovers the underlying shear fields and all of
the power spectra, including at levels well below the shot noise.
