Defining five-parameter Pearson distribution by moments

Question

How do you derive (or where can I find a derivation of) the equations to match the five parameters of a general Pearson distribution to a mean, variance, skewness, kurtosis, and condition that the function integrate to one?

I have been using the method described in the Mathematica documentation under "Examples" and then "Applications". It provides the equations but no derivation. The four equations for the moments are set up as:

eq[r_] := r*Subscript[b, 0]*Moment[r - 1] + (r + 1) Subscript[b, 1]*Moment[r] 
+ (r + 2) Subscript[b, 2]*Moment[r + 1] - Subscript[a, 1]*Moment[r + 1] 
- Subscript[a, 0]*Moment[r]

and

meq = Table[MomentConvert[eq[r], CentralMoment], {r, 0, 3}] /. 
    {Moment[1] -> μ, CentralMoment[2] -> σ^2, CentralMoment[3] -> Sqrt[Subscript[β, 1]] σ^3, 
        CentralMoment[4] -> Subscript[β, 2] σ^4}

These set up four linear equations in five unknowns.

I am interested in the equation added to fix the normalizing coefficient (so the distribution integrates to one):

meq = Join[meq, {Subscript[a, 0] + (12 μ Subscript[β, 1] + 
   2 μ (9 - 5 Subscript[β, 2]) + σ Sqrt[Subscript[β, 1]] (3 + Subscript[β, 2]))}];

The equations for the moments themselves are pretty easy to derive, but I have not figured out the fifth equation for normalization.

Answer 1

We solve this in Chapter 5 (Section 5.2 E: Higher Order Pearson-Style Families) of Mathematical Statistics with Mathematica. To the best of my knowledge, the 5 parameter Pearson had not been solved prior to the first edition of our book (if someone is aware of an earlier solution, I'd be pleased to hear of it). In particular, consider a Pearson-style system based upon a cubic polynomial. This will be the family of solutions $p(x)$ to the differential equation: $$\frac{dp(x)}{dx}=-\frac{(a+x) p(x)}{c_3 x^3+c_2 x^2+c_1 x+c_0}$$

After some work ... too long to show here, but related to boundary conditions etc ... one ends up with a recurrence relation in raw moments:

$$c_0 (-r) \acute{\mu }_{r-1}-c_2 (r+2) \acute{\mu }_{r+1}-c_3 (r+3) \acute{\mu }_{r+2}-c_1 (r+1) \acute{\mu }_r=-a \acute{\mu }_r-\acute{\mu }_{r+1}$$

Given the boundary condition $\acute{\mu }_0=1$, we enter this recurrence relation into Mathematica as:

Our objective is to find $a$, $c_0$, $c_1$, $c_2$ and $c_3$ in terms of$\acute{\mu }_r$. Putting $r = 0, 1, 2, 3, 4$ yields the required 5 equations (for the 5 unknowns) which we now solve simultaneously to yield the solution:

Z1 = Solve[Table[eqn2[r], {r, 0, 4}], {a, c0, c1, c2, c3}]

The solution is rather long, so I won't print the output here. However, if we work about the mean, taking $\acute{\mu }_1=0$, and $\acute{\mu }_r= \mu_r$ for $r>=2$, the solution reduces to:

which is comparatively compact. The text also provides examples with actual fitting of the third-order (cubic) Pearson, and comparison to the automated second-order (quadratic) Pearson family.