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

[josef.pktd at gmail.com]

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
-------------------------------------------