Abstract We detail a new method for performing robust Bayesian estimation of the three-dimensional spatial power spectrum of the Epoch of Reionization (EoR) from interferometric observations. The versatility of this technique allows us to present two approaches. First, when the observations span only a small number of independent spatial frequencies (k-modes), we sample directly from the spherical power spectrum coefficients that describe the EoR signal realization. Secondly, when the number of k-modes to be included in the model becomes large, we sample from the joint probability density of the spherical power spectrum and the signal coefficients, using Hamiltonian Monte Carlo methods to explore this high-dimensional (∼20 000) space efficiently. This approach has been successfully applied to simulated observations that include astrophysically realistic foregrounds in a companion publication (Sims et al. 2016). Here, we focus on explaining the methodology in detail and use simple foreground models to both demonstrate its efficacy and highlight salient features. In particular, we show that including an arbitrary flat spectrum continuum foreground that is 108 times greater in power than the EoR signal has no detectable impact on our parameter estimates of the EoR power spectrum recovered from the data.