  We propose a rigorous, non-perturbative, Bayesian framework which enables one
jointly to test Gaussianity and estimate the power spectrum of CMB
anisotropies. It makes use of the Hilbert space of an harmonic oscillator to
set up an exact likelihood function, dependent on the power spectrum and on a
set of parameters $\alpha_i$, which are zero for Gaussian processes. The latter
can be expressed as series of cumulants; indeed they perturbatively reduce to
cumulants. However they have the advantage that their variation is essentially
unconstrained. Any truncation(i.e.: finite set of $\alpha_i$) therefore still
produces a proper distribution - something which cannot be said of the only
other such tool on offer, the Edgeworth expansion. We apply our method to Very
Small Array (VSA) simulations based on signal Gaussianity, showing that our
algorithm is indeed not biased.
