  We introduce a new method for performing robust Bayesian estimation of the
three-dimensional spatial power spectrum at 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
realisation. Second, 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 ($\sim$ 20000) 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 $10^8$ 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.
