In this article, an efficient mathematical procedure is proposed for estimating a multivariate probability density function (mpdf). The method is based on statistical moments but does not require all samples. This is achieved by relating the probability distribution with these statistical moments. The novelty of the method is the calculation of coefficients associated with the mpdf on a suitable Riesz basis chosen for (Formula presented.). The chosen base has several advantages. First, once the coefficients are known to calculate the mpdf, it does not require any matrix inversion. Second, these coefficients can be precalculated and stored. Third, the method is modular, which makes it possible to parallelize the process. These advantages lead us to a very accurate and economical procedure from a cost perspective of computational efforts and time savings. Different numerical examples are presented to prove the effectiveness of our proposed methodology and are compared with more traditional approaches, such as the kernel density method.