  We consider the shape of the posterior distribution to be used when fitting
cosmological models to power spectra measured from galaxy surveys. At very
large scales, Gaussian posterior distributions in the power do not approximate
the posterior distribution $\mathcal{P}_R$ we expect for a Gaussian density
field $\delta_k$, even if we vary the covariance matrix according to the model
to be tested. We compare alternative posterior distributions with
$\mathcal{P}_R$, both mode-by-mode and in terms of expected measurements of
primordial non-Gaussianity parameterised by $f_\mathrm{NL}$. Marginalising over
a Gaussian posterior distribution $\mathcal{P}_f$ with fixed covariance matrix
yields a posterior mean value of $f_\mathrm{NL}$ which, for a data set with the
characteristics of Euclid, will be underestimated by $\triangle
f_\mathrm{NL}=0.4$, while for the data release 9 (DR9) of the Sloan Digital Sky
Survey (SDSS)-III Baryon Oscillation Spectroscopic Survey (BOSS) it will be
underestimated by $\triangle f_\mathrm{NL}=19.1$. Adopting a different form of
the posterior function means that we do not necessarily require a different
covariance matrix for each model to be tested: this dependence is absorbed into
the functional form of the posterior. Thus, the computational burden of
analysis is significantly reduced.
