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.

• The first $k$ moments define the first $k$ derivatives of the characteristic function at zero: $E[X^k] = (-i)^k\phi_X^{(k)}(0)$. So you know the first $k$ terms of the characteristic function's Taylor expansion around zero. You may then be able to use the inversion theorems to derive the density. – Stephan Kolassa Oct 3 '15 at 14:07
• Thanks @StephanKolassa - any chance for an extended answer / an R code example? – Tal Galili Oct 3 '15 at 15:37
• en.wikipedia.org/wiki/Maximum_entropy_probability_distribution suggests a general method. – whuber Oct 3 '15 at 17:40
• Dear @whuber, could you please suggest an R code example? (also, does this go with wolfies's answer?) – Tal Galili Oct 3 '15 at 21:27
• This is a completely different approach from that answer. – whuber Oct 3 '15 at 22:21

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.

• Dear @wolfies - thank you for your answer! If I understand you correctly, the Gram-Charlier expansions is more in line with what I am looking for (although the more generalized Pearson distribution is interesting to know about). I looked at your book (chapter 5, starting from page 175), and see you indeed give there a detailed description (with also mentions of how to deal with estimated moments, which is my case). The only thing is that I can't use your code (since I am an R user). Thanks for your answer (and also for your book which seems impressive and interesting in general) – Tal Galili Oct 3 '15 at 20:43
• Just found an R package to deal with the various methods: cran.us.r-project.org/web/packages/PDQutils/vignettes/… – Tal Galili Oct 3 '15 at 20:48