[SciPy-User] rv_frozen when using gamma function
Bruno Santos
bacmsantos at gmail.com
Mon Mar 19 12:34:48 EDT 2012
[quote]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
[/quote]
Guys you are starting to overwhelm my knowledge :). But you are correct the
only way I can solve it by assuming that d is given and fixed. It is not
ideally but unfortunately I haven't found a way to get it experimentally so
will have to pretend I know and use several values.
And Skipper thank you very much some times the solution is so obvious I
have problems in seeing it. The -1 did the trick.
Thank you very much for all the help. I am really learning a lot from this
thread.
Bruno
On 19 March 2012 16:28, <josef.pktd at gmail.com> wrote:
>
>
> On Mon, Mar 19, 2012 at 12:23 PM, Skipper Seabold <jsseabold at gmail.com>wrote:
>
>> 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.
>>
>
> 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 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
>>>
>>>
>>
>> _______________________________________________
>> 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/8e322c9f/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/8e322c9f/attachment.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/8e322c9f/attachment-0001.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/8e322c9f/attachment-0002.png>
More information about the SciPy-User
mailing list