Search code examples
pythonmatlabscipy

Gammainc difference in Python and Matlab

I am transfering code between Matlab and Python. However when I use scipy.special.gammainc I get different values than when I use gammainc in Matlab.

My Python code is:

from scipy.special import gammainc                                        
gammainc([1, 2, 3, 4], 5)                                                
array([0.99326205, 0.95957232, 0.87534798, 0.73497408])

However in Matlab I have:

gammainc([1,2,3,4],5,'lower')                                      
[0.00365984682734371, 0.0526530173437111, 0.184736755476228, 0.371163064820127]

I am trying to change my Python version to give me the same version as the Matlab code. I have tried instead:

from scipy.special import gammaincc                                   
gammaincc([1, 2, 3, 4], 5)                                               
array([0.00673795, 0.04042768, 0.12465202, 0.26502592])

This is closer but it is still not what I am looking for.

Solution

  • There are two issues with your code:

    • Matlab's gammainc function expects the input arguments in reverse order compared to the usual mathematical notation. So it expects the integral limit first, then the exponent.
    • Scipy's gammaincc is the upper incomplete Gamma function, not the lower. For the lower version you need to use gammainc.

    So, to reproduce Matlab's result in Python with Scipy, use

    from scipy.special import gammainc
print(gammainc(5, [1, 2, 3, 4]))

    This gives

    [0.00365985 0.05265302 0.18473676 0.37116306]