  We derive and implement a full Bayesian large scale structure inference
method aiming at precision recovery of the cosmological power spectrum from
galaxy redshift surveys. Our approach improves over previous Bayesian methods
by performing a joint inference of the three dimensional density field, the
cosmological power spectrum, luminosity dependent galaxy biases and
corresponding normalizations. We account for all joint and correlated
uncertainties between all inferred quantities. Classes of galaxies with
different biases are treated as separate sub samples. The method therefore also
allows the combined analysis of more than one galaxy survey.
  In particular, it solves the problem of inferring the power spectrum from
galaxy surveys with non-trivial survey geometries by exploring the joint
posterior distribution with efficient implementations of multiple block Markov
chain and Hybrid Monte Carlo methods. Our Markov sampler achieves high
statistical efficiency in low signal to noise regimes by using a deterministic
reversible jump algorithm. We test our method on an artificial mock galaxy
survey, emulating characteristic features of the Sloan Digital Sky Survey data
release 7, such as its survey geometry and luminosity dependent biases. These
tests demonstrate the numerical feasibility of our large scale Bayesian
inference frame work when the parameter space has millions of dimensions.
  The method reveals and correctly treats the anti-correlation between bias
amplitudes and power spectrum, which are not taken into account in current
approaches to power spectrum estimation, a 20 percent effect across large
ranges in k-space. In addition, the method results in constrained realizations
of density fields obtained without assuming the power spectrum or bias
parameters in advance.
