[SciPy-User] rv_frozen when using gamma function

Mon Mar 19 12:14:58 EDT 2012

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?

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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20120319/9f7e25d4/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/9f7e25d4/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/9f7e25d4/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/9f7e25d4/attachment-0002.png>