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.