  In this work we investigate the multivariate statistical description of the
matter distribution in the nonlinear regime. We introduce the multivariate
Edgeworth expansion of the lognormal distribution to model the cosmological
matter field. Such a technique could be useful to generate and reconstruct
three-dimensional nonlinear cosmological density fields with the information of
higher order correlation functions. We explicitly calculate the expansion up to
third order in perturbation theory making use of the multivariate Hermite
polynomials up to sixth order. The probability distribution function for the
matter field includes at this level the two-point, the three-point and the
four-point correlation functions. We use the hierarchical model to formulate
the higher order correlation functions based on combinations of the two-point
correlation function. This permits us to find compact expressions for the
skewness and kurtosis terms of the expanded lognormal field which can be
efficiently computed. The method is, however, flexible to incorporate arbitrary
higher order correlation functions which have analytical expressions. The
applications of such a technique can be especially useful to perform
weak-lensing or neutral hydrogen 21 cm line tomography, as well as to directly
use the galaxy distribution or the Lyman-alpha forest to study structure
formation.
