[SciPy-User] stats.distributions moments and expect - another round

[josef.pktd at gmail.com]

On Sat, Dec 4, 2010 at 9:27 AM,  <josef.pktd at gmail.com> wrote:
> I spend most of a day fixing the expect function for stats
> distributions and checking mean, variance, skew and kurtosis, and
> trying to get it to work for almost all distributions.
>
> expect which uses integrate quad doesn't always work well, but the
> results look mostly good except for some fat-tailed distributions.
>
> Below are the comparisons between the stats method of the
> distributions and the outcome of expect. The first failure for mean is
> at 4 decimals. The tables are discrepancies at two decimals:
>
> invgamma might still have some numerical integration problems, that I
> haven't figured out yet. Some other ones still look like incorrect
> formulas for skew and kurtosis.
>
> As an aside: I think, finally, that we can add the doc templates to
> many distributions, so we have a way of pointing out differences in
> parameterization and known numerical problems. For example, I'm still
> not sure whether dist.stats contains the correct information about
> whether the moments don't exist or are infinite.
>
> Josef
>
>
> Variance
>
>>>> print SimpleTable(var_, headers=['distname', 'diststats', 'expect'])
> ==========================================
>  distname    diststats         expect
> ------------------------------------------
>   pareto   1.60340572053   1.57246536012
> tukeylambda 0.304764722791 0.0268724542194
> fatiguelife   884942.25      884940.75504
>     t      3.69051720817   3.66604272097
>  powerlaw  0.858574961358 0.0641248702779
>  invgamma  13.1319516251    5.5288716506
>   rdist    1.78002848501   0.526315790145
> ------------------------------------------
>
> Skew
>
>>>> print SimpleTable(skew, headers=['distname', 'diststats', 'expect'])
> ===========================================
>  distname     diststats         expect
> -------------------------------------------
>   mielke    7.59540085257   6.91088104518
>    fisk     38.7938857832   12.9246301568
>  foldnorm   0.971407236222  0.202188769695
>  gilbrat    6.18487713863   6.17333365849
>  loglaplace  16.9237038681   11.2535633383
> fatiguelife 0.0408718910942  3.93239048725
>  powerlaw  -0.906960466124 -0.420181302329
>    ncf      40747519832.7   8.94481163856
>     f       1.93130205529   1.80641432186
>  invgamma  -0.477729689377  655.942282864
> -------------------------------------------
>
> Kurtosis
>
>>>> print SimpleTable(kurt, headers=['distname', 'diststats', 'expect'])
> ===========================================
>  distname     diststats         expect
> -------------------------------------------
>   mielke    -149.405089743  362.442518204
>    fisk     -224.659270348  2099.30241789
>  foldnorm   2.70517285483  -0.294828589688
> tukeylambda  -2.98365209914 -0.897302898918
>  dweibull   1.90893020344   -1.06484211833
>  gilbrat    110.936392176   107.859563214
>  loglaplace  -164.332555303   1330.2395564
>  genpareto   14.8285714286   14.8119563403
>  lognorm    81.1353811489   79.6180870122
>    burr     112616.270172   6.21265707889
>    ncf     -239984516633.0  13492.9902378
>     f        7.9138697318   7.06539159862
>    nct      -409040.407062  0.605963897342
>  invgamma    -2.866573514   2116889.58176
>   rdist     -2.56785479799  -1.53846154515
> -------------------------------------------
>

one down: invgamma is correct if a>4, requirement for kurtosis to exist

Josef