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

josef.pktd@gmai... josef.pktd@gmai...
Sat Dec 4 10:29:15 CST 2010


On Sat, Dec 4, 2010 at 10:35 AM,  <josef.pktd@gmail.com> wrote:
> On Sat, Dec 4, 2010 at 9:27 AM,  <josef.pktd@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

I'm giving up, staring at the kurtosis of nct and burr without the
direct references looks like more fun than I have time for. I think I
have a patch for ncf.

Any volunteers, or references ?
At least we are down to a reasonably short list that needs checking or
bugfixing.

Josef

>
> Josef
>


More information about the SciPy-User mailing list