How to fit an approximate PDF (i.e.: density estimation) using the first k (empirical) moments? I have a situation where I am able to estimate (the first) $k$ moments of a data-set, and would like to use it to produce an estimation of the density function.
I already came across the Pearson distribution, but realized it relies only on the first 4 moments (with some restrictions on the possible combinations of moments).
I also understand that any finite set of moments is not enough to "pin down" a specific distribution, when not using more assumptions. However, I would still like for a more general class of distributions (other than the Pearson family of distributions). Looking at other questions, I couldn't find such a distribution (see: here, here, here, here, here, and here).
Is there some ("simple") generalized family of distribution that can be defined for any set of $k$ moments? (maybe a set of transformations that can take a standard normal distribution and transforms it until it confirms with all set of $k$ moments)
(I don't care much if we assume the other $k+1\ldots\infty$ moments are 0 or not)
Thanks.
p.s.: I would be happy for an extended example. Preferably with an R code example.
 A: Method 1: Higher-order Pearson systems
The Pearson system is, by convention, taken to be the family of solutions $p(x)$ to the differential equation:
$$
\frac{d p (x)}{dx} \; = \; -\frac{(a+x) }{c_0 + c_1 x + c_2 x^2} \; p(x)
$$
where the four Pearson parameters $(a, c_0, c_1, c_2)$ can be expressed in terms of the first four moments of the population. 
Instead of basing the Pearson system upon the quadratic $c_0 + c_1 x + c_2 x^2$, one can instead consider using higher order polynomials as the foundation stone. So, for example, one can consider a Pearson-style system based upon a cubic polynomial. This will be the family of solutions $p(x)$ to the differential equation:
$$
\frac{d p(x)}{dx} \; = \; -\frac{(a+x) }{c_0 + c_1 x + c_2 x^2 + c_3 x^3} \; p (x)
$$
which yields the solution:

I solved this for fun some time back (having the same thought train as the OP): the derivation and solution is given in Chapter 5 of our book; if interested, a free download is available here:
http://www.mathstatica.com/book/bookcontents.html
Note that whereas the second-order (quadratic) Pearson family can be expressed in terms of the first 4 moments, the third-order (cubic) Pearson-style family requires the first 6 moments.
Method 2: Gram-Charlier expansions
Gram-Charlier expansions are also discussed in the same Chapter 5 (see section 5.4) ... and also allow one to construct a fitted density, based on arbitrarily large $k^{th}$ moments. As the OP suggests, the Gram-Charlier expansion expresses the fitted pdf as a function of a series of derivatives of the standard Normal pdf, known as Hermite polynomials. The Gram-Charlier coefficients are solved as a function of the population moments ... and the bigger the expansion, the more moments required. You may also wish to look at related Edgeworth expansions.
Population moments or sample moments??
For the Pearson-style system:  If the moments of the population are known, then using higher moments should unambiguously yield a better fit. If, however, the observed data is a random sample drawn from the population, there is a trade-off: a higher order polynomial implies that higher order moments are required, and the estimates of the latter may be unreliable (have high variance), unless the sample size is 'large'. In other words, given sample data, fitting using higher moments can become 'unstable' and produce inferior results. The same is true for Gram-Charlier expansions: adding an extra term can actually yield a worse fit, so some care is required.
