Pearson Distribution moment derivation

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 (PearsonDistribution—Wolfram 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] -> \[Mu], CentralMoment[2] -> \[Sigma]^2, CentralMoment[3] -> Sqrt[Subscript[\[Beta], 1]] \[Sigma]^3, CentralMoment[4] -> Subscript[\[Beta], 2] \[Sigma]^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 \[Mu] Subscript[\, 1] + 2 \[Mu] (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.

Sorry for answering this so late … I have only just seen the question (by chance).
The Pearson family is conventionally defined in terms of 4 parameters (not 5). Solving the Pearson coefficients in terms of moments involves assuming the density tends to zero at the extremum of the tails, solving a recurrence relation, and is a little bit of work, but Mathematica does it very nicely indeed. See section 5.2D of our book, “Mathematical Statistics with Mathematica” which provides the derivation; if interested, a free download is available here:

Section 5.2E of Chapter 5 extends to a 5-parameter model based on a third-order (cubic) polynomial, instead of the standard second-order (quadratic) model … and the same method works just as nicely.