[SciPy-dev] Problem with F distribution, or with me?

josef.pktd@gmai... josef.pktd@gmai...
Wed Aug 13 13:44:14 CDT 2008


I think, stats.loggamma.rvs is wrong or uses a definition that I cannot
figure out.
It was still bugging me, that I did not figure out what is going on with the
loggamma random variables.

When I compare it with the log transformation of a gamma random variable and
with the loggamma distribution in R, then it looks llike as if the
stats.loggamma.rvs were only correct for a parameter of 2.

In "R", I get the same mean and variance as for
np.log(stats.gamma.rvs(k,size=1000000))
and neither R nor np.log(stats.gamma.rvs(...)) have the same domain
restriction as stats.loggamma.rvs.

So, to me it seems that there is something fishy with stats.loggamma.rvs.

used package in R: 'VGAM'

Below are some comparison of the loggamma random variables generated by
 * stats.loggamma.rvs
 * np.log(stats.gamma.rvs(...))
 * VGAM.rlgamma in "R"

Josef

Kolmogorov tests
------------------------

for parameter = 2 it looks ok

>>> stats.kstest(np.log(stats.gamma.rvs(2,size=10000)),'loggamma',args=(2,))
(0.010148906122060208, array(0.12659359569633333))
>>>
c=2;stats.kstest(np.log(stats.gamma.rvs(c,size=10000)),'loggamma',args=(c,))
(0.0058379222756740345, array(0.50383387183308759))

strange results with c different from 2

>>>
c=2.5;stats.kstest(np.log(stats.gamma.rvs(c,size=10000)),'loggamma',args=(c,))
(0.24775311834187796, array(0.0))
>>>
c=2.5;stats.kstest(stats.loggamma.rvs(c,size=10000),'loggamma',args=(c,))
(1.0, array(0.0))
>>>
c=1.5;stats.kstest(stats.loggamma.rvs(c,size=10000),'loggamma',args=(c,))
(0.0069725721519633965, array(0.3764490603546865))
>>>
c=1.5;stats.kstest(np.log(stats.gamma.rvs(c,size=10000)),'loggamma',args=(c,))
(0.12849906523207333, array(0.0))


Comparing mean and variance with "R"
-------------------------------------------------------

for k=2: is ok

>>> np.mean(stats.loggamma.rvs(2,size=10000))
0.42482669940021239
>>> np.mean(np.log(stats.gamma.rvs(2,size=100000)))
0.42284709204863447
in R
> mean(rlgamma(100000, location=0, scale=1, k=2))
[1] 0.4224025


>>> np.var(np.log(stats.gamma.rvs(2,size=10000)))
0.6449444777702128
>>> np.var(stats.loggamma.rvs(2,size=100000))
0.64758634320780184
in R:
> var(rlgamma(100000, location=0, scale=1, k=2))
[1] 0.6360736


for k=5:

>>> np.mean(stats.loggamma.rvs(5,size=1000000))
1.#QNAN
>>> np.mean(np.log(stats.gamma.rvs(5,size=1000000)))
1.506037882165763
in R:
> mean(rlgamma(1000000, location=0, scale=1, k=5))
[1] 1.506325

>>> np.var(stats.loggamma.rvs(5,size=1000000))
1.#QNAN
>>> np.var(np.log(stats.gamma.rvs(5,size=1000000)))
0.22151419947589021
> var(rlgamma(1000000, location=0, scale=1, k=5))
[1] 0.2212415

for other k:

>>> np.var(stats.loggamma.rvs(1.5,size=1000000))
0.78615045130829264
>>> np.var(np.log(stats.gamma.rvs(1.5,size=1000000)))
0.93399967614023915
in R:
> var(rlgamma(1000000, location=0, scale=1, k=1.5))
[1] 0.9332276


>>> np.var(np.log(stats.gamma.rvs(10,size=1000000)))
0.10528311964943039
in R:
> var(rlgamma(1000000, location=0, scale=1, k=10))
[1] 0.1052298
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-dev/attachments/20080813/02e79615/attachment.html 


More information about the Scipy-dev mailing list