Power spectrum estimation from high-resolution maps by Gibbs sampling

We revisit a recently introduced power spectrum estimation technique based on Gibbs sampling, with the goal of applying it to the high-resolution WMAP data. In order to facilitate this analysis, a number of sophistications have to be introduced, each of which is discussed in detail. We have implemented two independent versions of the algorithm to cross-check the computer codes, and to verify that a particular solution to any given problem does not affect the scientific results. We then apply these programs to simulated data with known properties at intermediate (N_side = 128) and high (N_side = 512) resolutions, to study effects such as incomplete sky coverage and white vs. correlated noise. From these simulations we also establish the Markov chain correlation length as a function of signal-to-noise ratio, and give a few comments on the properties of the correlation matrices involved. Parallelization issues are also discussed, with emphasis on real-world limitations imposed by current super-computer facilities. The scientific results from the analysis of the first-year WMAP data are presented in a companion letter.