[SciPy-user] bug in stats.pdfapprox/stats.pdf_moment and new Gram-Charlier distribution

[josef...@gmail.com]

I wanted to create a new distribution by wrapping stats.pdfapprox and
stats.pdf_moment in stats.morestats.py.

However, these two function do not create a normal expansion if the
first four moments are given. The inner loop in stats.pdf_moment is
never entered when four moments are given. As a consequence the pdf
that is returned is the unexpanded normal distribution. (There is also
a small mistake in stats.pdfapprox in how the variance is calculated.)

I didn't find any information which type of expansion is used by
pdf_moment. I assume it is Gram-Charlier, but I didn't find any
formulas to make sense out of the inner loop that calculates the
coefficients (for multiplying with the Hermite polynomials).

If someone could provide a (understandable) reference for these
calculations or figure out what the loop is supposed to do, then we
could correct the expansion for the general case.

Since I couldn't fix `pdf_moment`, I wrote a new function that
calculates the pdf for the Gram-Charlier expansion when the first four
moments (or mean, variance, skew, kurtosis) are given. This uses the
explicit formula for this expansion, and doesn't allow for higher
order expansion.

pdf_mvsk: get pdf of G-Ch normal expansion using mean, variance, skew,
and excess kurtosis

This I wrapped in a subclass of _distributionsrv_continuous: NormExpan_gen
It works in the examples that I tried but is not fully tested or cleaned up yet.

attachment:
* try_pdfapprox.py shows problem with current function
* distr_gch.py new expansion pdf, and NormExpan distribution

I also wrote a skew normal and skew t distribution (as defined by
Azzalini, A. & Capitanio, A., univariate only), which is not attached.

Josef