Using scipy.integrate.quad to integrate gamma function

Question

How do I find the integral of Gamma(x) = x^(a-1) * e^-x from 0 to infinity when a = 3/2 using scipy.integrate.quad?

So far I've tried some code but I am being told that "quad: The first argument is not callable".

This is the code I've used:

import numpy as np 
import math as m 
import matplotlib.pyplot as plt
from math import exp
from scipy import *
import scipy.integrate

def f(a):
    return x**(a-1)*exp(-x)

a = 3/2

print scipy.integrate.quad(f(a), 0, inf)

Answer 1

You could define f to take an argument for a as follows:

def f(x, a):
    return x**(a-1) * np.exp(-x)

Note that it's necessary for f to have the first argument as x since this is the variable we'll be integrating with respect to.

To integrate, you can then pass the value of a to f using args:

>>> a = 3./2
>>> scipy.integrate.quad(f, 0, np.inf, args=a)
(0.8862269254536111, 9.077554263825505e-10)

A couple of other points:

• Your original error was due to the fact that quad needs a function (or callable) object as its first argument. Writing f(...) passes the value returned from the function, not the function object itself.
• In Python 2, the expression 3/2 will give the integer 1, not the float 1.5. You need to make one of the numbers a float to trigger true division (Python 3 does not have this problem).