# [SciPy-User] rv_frozen when using gamma function

josef.pktd@gmai... josef.pktd@gmai...
Mon Mar 19 11:28:14 CDT 2012

```On Mon, Mar 19, 2012 at 12:23 PM, Skipper Seabold <jsseabold@gmail.com>wrote:

> On Mon, Mar 19, 2012 at 12:14 PM, Bruno Santos <bacmsantos@gmail.com>wrote:
>
>> I believe the formula I have is accurate I checked it personally and also
>> have it checked by two mathematicians in the lab and they come up with the
>> same results. I left my notebook where I performed the transformations home
>> so don't completely remember but I believe you can simply things to get rid
>> of some of the parameters.
>>
> dicerAcc is a scalar as you mentioned.
>> I managed to implement the function in python now and it is giving the
>> same results as in R my question how to maximize it still remains though.
>> Is it possibly to maximize a function rather than minimize it in Python?
>>
>
> Ok, then I guess my math is faulty. I only looked quickly and don't see
> the other close parens in the formula.
>

I didn't check the parens, but to me it just like the negative binomial but
I think only the ratio

p = mu/(beta+mu)
1-p = beta/(beta+mu)

are identified, unless a parameter is held fixed.

negative binomial is a poisson-gamma mixture, but I only found the p, 1-p
parameterization

Josef

>
> To maximize put a negative in front of the function.
>
>
>
>>
>>
>> On 15 March 2012 15:21, Skipper Seabold <jsseabold@gmail.com> wrote:
>>
>>> On Thu, Mar 15, 2012 at 11:07 AM, Bruno Santos <bacmsantos@gmail.com>wrote:
>>>
>>>> Thank you all very much for the replies that was exactly what I wanted.
>>>> I am basically trying to get the parameters for a
>>>> gamma-poisson distribution. I have the R code from a
>>>> previous collaborator just trying to write a native function in python
>>>> rather than using the R code or port it using rpy2.
>>>
>>>
>>> Oh, fun.
>>>
>>>
>>>> The function is the following:
>>>> [image: Inline images 1]
>>>> where f(b,d) is a function that gives me a probability of a certain
>>>> position in the vector to be occupied and it depends on b (the position)
>>>> and d (the likelihood of making an error).
>>>> So the likelihood after a few transformations become:
>>>>
>>>> [image: Inline images 2]
>>>> Which I then use the loglikelihood and try to maximise it using an
>>>> optimization algorithm.
>>>> [image: Inline images 3]
>>>> The R code is as following:
>>>> alphabeta<-function(alphabeta,x,dicerAcc)
>>>> {
>>>>   alpha <-alphabeta[1]
>>>>   beta <-alphabeta[2]
>>>>   if (any(alphabeta<0))
>>>>     return(NA)
>>>>   sum((alpha*log(beta) + lgamma(alpha + x) + x * log(dicerAcc) -
>>>> lgamma(alpha) - (alpha + x) * log(beta+dicerAcc) - lfactorial(x))[dicerAcc
>>>> > noiseT])
>>>>
>>>
>>> From a quick (distracted) look (so I could be wrong)
>>>
>>> Should this be alpha^2*log(beta) ? +lgamma(alpha) ? And lfactorial(x)
>>> should still be +lgamma(alpha)*lfactorial(x) ? And dicerAcc a scalar
>>> integer I take it?
>>>
>>>
>>>>
>>>> #sum((alpha*log(beta)+(lgamma(alpha+x)+log(dicerError^x))-(lgamma(alpha)+log((beta+dicerError)^(alpha+x))+lfactorial(x)))[dicerError
>>>> != 0])
>>>> }
>>>> x and dicerAcc are known so the I use the optim function in R
>>>> ab <- optim(c(1,100), alphabeta, control=list(fnscale=-1), x = x,
>>>> dicerAcc = dicerAcc)\$par
>>>>
>>>> Is there any equivalent function in Scipy to the optim one?
>>>>
>>>> On 14 March 2012 17:05, Bruno Santos <bacmsantos@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)))
>>>>>
>>>>> 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@scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>
>>>>
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/7e0c90fd/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 4401 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/7e0c90fd/attachment-0003.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 1620 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/7e0c90fd/attachment-0004.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: image/png
Size: 3193 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/7e0c90fd/attachment-0005.png
```

More information about the SciPy-User mailing list