[SciPy-User] rv_frozen when using gamma function

Skipper Seabold jsseabold@gmail....
Thu Mar 15 10:21:51 CDT 2012


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


More information about the SciPy-User mailing list