[SciPy-User] rv_frozen when using gamma function

[josef.pktd at gmail.com]

On Wed, Mar 14, 2012 at 1:05 PM, Bruno Santos <bacmsantos at gmail.com> wrote:
> I am trying to write a script to do some maximum likelihood parameter
> estimation of a function. But when I try to use the gamma function I get:
> gamma(5)
> Out[5]: <scipy.stats.distributions.rv_frozen at 0x7213710>
>
> I thought it might have been a problem solved already on the new
> distribution but even after installing the last scipy version I get the same
> problem.
> The test() after installation is also failing with the following
> information:
> Running unit tests for scipy
> NumPy version 1.5.1
> NumPy is installed in /usr/lib/pymodules/python2.7/numpy
> SciPy version 0.10.1
> SciPy is installed in /usr/local/lib/python2.7/dist-packages/scipy
> Python version 2.7.2+ (default, Oct  4 2011, 20:06:09) [GCC 4.6.1]
> nose version 1.1.2
> ...
> ...
> ...
> AssertionError:
> Arrays are not almost equal
>  ACTUAL: 0.0
>  DESIRED: 0.5
>
> ======================================================================
> FAIL: Regression test for #651: better handling of badly conditioned
> ----------------------------------------------------------------------
> Traceback (most recent call last):
>   File
> "/usr/local/lib/python2.7/dist-packages/scipy/signal/tests/test_filter_design.py",
> line 34, in test_bad_filter
>     assert_raises(BadCoefficients, tf2zpk, [1e-15], [1.0, 1.0])
>   File "/usr/lib/pymodules/python2.7/numpy/testing/utils.py", line 982, in
> assert_raises
>     return nose.tools.assert_raises(*args,**kwargs)
> AssertionError: BadCoefficients not raised
>
> ----------------------------------------------------------------------
> Ran 5103 tests in 47.795s
>
> FAILED (KNOWNFAIL=13, SKIP=28, failures=3)
> Out[7]: <nose.result.TextTestResult run=5103 errors=0 failures=3>
>
>
> My code is as follows:
> from numpy import array,log,sum,nan
> from scipy.stats import gamma
> from scipy import factorial, optimize
>
> #rinterface.initr()
> #IntSexpVector = rinterface.IntSexpVector
> #lgamma = rinterface.globalenv.get("lgamma")
>
> #Implementation for the Zero-inflated Negative Binomial function
> def alphabeta(params,x,dicerAcc):
>     alpha = array(params[0])
>     beta = array(params[1])
>     if alpha<0 or beta<0:return nan
>     return sum((alpha*log(beta)) + log(gamma(alpha+x)) + x * log(dicerAcc) -
> log(gamma(alpha)) - (alpha+x) * log(beta+dicerAcc) - log(factorial(x)))

I guess what you want her is scipy.special.gamma which is the gamma
function, not the gamma distribution

loglikelihood of negative binomial is also in statsmodels.discrete if
you want to compare notes

Josef


>
> if __name__=='__main__':
>     x =
> array([123,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,104,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,24,1,0,0,0,0,0,0,0,2,0,0,4,0,0,0,0,0,0,0,0,12,0,0])
>     dicerAcc = array([1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
> 0.048750000000000002,0.90085000000000004, 0.0504, 0.0, 0.0, 0.0, 0.0, 0.0,
> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0023,
> 0.089149999999999993, 0.81464999999999999, 0.091550000000000006,
> 0.0023500000000000001, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
> 0.0, 0.0, 0.0, 0.0, 0.0, 0.00020000000000000001, 0.0061000000000000004,
> 0.12085, 0.7429, 0.12325, 0.0067000000000000002, 0.0, 0.0, 0.0, 0.0, 0.0,
> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.00020000000000000001,
> 0.012500000000000001, 0.14255000000000001, 0.68159999999999998,
> 0.14979999999999999, 0.012999999999999999])
>     optimize.()
>
>
> Am I doing something wrong or is this a known problem?
>
> Best,
> Bruno
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>