[SciPy-User] rv_frozen when using gamma function

Mon Mar 19 12:23:20 EDT 2012

On Mon, Mar 19, 2012 at 12:14 PM, Bruno Santos <bacmsantos at 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.

To maximize put a negative in front of the function.

>
>
> On 15 March 2012 15:21, Skipper Seabold <jsseabold at gmail.com> wrote:
>
>> On Thu, Mar 15, 2012 at 11:07 AM, Bruno Santos <bacmsantos at 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 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)))
>>>>
>>>> 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
>>>
>>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at 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/adafd28e/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 4401 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/adafd28e/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 1620 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/adafd28e/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 3193 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/adafd28e/attachment-0002.png>