  Non-linear gravitational collapse introduces non-Gaussian statistics into the
matter fields of the late Universe. As the large-scale structure is the target
of current and future observational campaigns, one would ideally like to have
the full probability density function of these non-Gaussian fields. The only
viable way we see to achieve this analytically, at least approximately and in
the near future, is via the Edgeworth expansion. We hence rederive this
expansion for Fourier modes of non-Gaussian fields and then continue by putting
it into a wider statistical context than previously done. We show that in its
original form, the Edgeworth expansion only works if the non-Gaussian signal is
averaged away. This is counterproductive, since we target the
parameter-dependent non-Gaussianities as a signal of interest. We hence alter
the analysis at the decisive step and now provide a roadmap towards a
controlled and unadulterated analysis of non-Gaussianities in structure
formation (with the Edgeworth expansion). Our central result is that, although
the Edgeworth expansion has pathological properties, these can be predicted and
avoided in a careful manner. We also show that, despite the non-Gaussianity
coupling all modes, the Edgeworth series may be applied to any desired subset
of modes, since this is equivalent (to the level of the approximation) to
marginalising over the exlcuded modes. In this first paper of a series, we
restrict ourselves to the sampling properties of the Edgeworth expansion,
i.e.~how faithfully it reproduces the distribution of non-Gaussian data. A
follow-up paper will detail its Bayesian use, when parameters are to be
inferred.
