[SciPy-User] rv_frozen when using gamma function

Bruno Santos bacmsantos@gmail....
Thu Mar 15 10:07:25 CDT 2012


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.
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])

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


More information about the SciPy-User mailing list