Defining five-parameter Pearson distribution by moments

Defining five-parameter Pearson distribution by moments

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
(http://reference.wolfram.com/mathematica/ref/PearsonDistribution.html)
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 for Mathematica.
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[\[Beta], 1] +
2 \[Mu] (9 - 5 Subscript[\[Beta], 2]) + \[Sigma] Sqrt[
Subscript[\[Beta], 1]] (3 + Subscript[\[Beta], 2]))}];
The equations for the moments themselves are pretty easy to derive, but I
have not figured out the fifth equation for normalization.