Plotting the gamma function in Python

The gamma function is defined for $x > 0$, $x \in \mathbb{R}$ as follows $\Gamma(x) = \int_{0}^{\infty} t^{x-1} e^{-t}\ dt$
This integral is improper at both its ends ($0$ and $\infty$) but it is convergent for every $x > 0$.
Also, when $x=1$ it could be read just as $\Gamma(1) = \int_{0}^{\infty} e^{-t}\ dt$

See also: Gamma Function - Wikipedia

In [1]:

import math as mt
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sp

In [2]:

# gm = np.vectorize(mt.gamma)
# fact = np.vectorize(mt.factorial)

gm = sp.gamma # the gamma function from scipy (already vectorized)
fact = sp.factorial # the factorial function from scipy (already vectorized)

In [3]:

X_MIN = 0.008
X_MAX = 6.1

x = np.linspace(X_MIN, X_MAX, 200)
y = gm(x)

k = np.arange(1, X_MAX)

fig, ax = plt.subplots(figsize=(10,6))

ax.plot(x, y, linewidth=2.0, label='$\Gamma(x)$ for $0.008 < x < 6.1$')
ax.plot(k, fact(k-1), 'ro', alpha=0.6,
         label='(k-1)! for k = 1, 2, 3, ...')

ax.set(xlim=(X_MIN-0.5, X_MAX+1))

plt.xlabel('X')
plt.ylabel('Y', rotation=0)

plt.legend(loc='lower right')
plt.show()

In [4]:

k

Out[4]:

array([1., 2., 3., 4., 5., 6.])

In [5]:

gm(0.001), gm(0.01), gm(0.1), gm(0.5), gm(1), gm(2), gm(3), gm(4), gm(5), gm(6)

Out[5]:

(999.4237724845956,
 99.43258511915059,
 9.513507698668732,
 1.7724538509055159,
 1.0,
 1.0,
 2.0,
 6.0,
 24.0,
 120.0)

In [6]:

# This demonstrates the important property that   
# gamma(1/2) is equal to sqrt(pi)   

abs(mt.sqrt(np.pi) - gm(0.5))

Out[6]:

0.0

In [7]:

X_MIN = 0.008
X_MAX = 1.01

x = np.linspace(X_MIN, X_MAX, 200)
y = gm(x)

k = np.arange(1, X_MAX)

fig, ax = plt.subplots(figsize=(10,6))

ax.plot(x, y, linewidth=2.0, label='$\Gamma(x)$ for $0.008 < x < 1.01$')

ax.set(xlim=(X_MIN-0.1, X_MAX+0.1))

plt.xlabel('X')
plt.ylabel('Y', rotation=0)

plt.legend(loc='upper right')
plt.show()

The lines of code below demonstrate several properties of the gamma function

$\Gamma\left(1/2\right) = \sqrt{\pi}$

Also, it is true that

$\Gamma\left((2m+1)/2\right) = \sqrt{\pi} \cdot \frac{1 \cdot 3 \dots (2m-1)}{2^m}$

for any natural number $m = 1,2,3, \dots$

Another important property is that

$\Gamma\left(n\right) = (n-1)!$

for any natural number $m = 1,2,3, \dots$

In [8]:

abs(gm(1.5) - mt.sqrt(np.pi) * 1 / 2) < 0.001

Out[8]:

True

In [9]:

abs(gm(2.5) - mt.sqrt(np.pi) * 3 / 4) < 0.001

Out[9]:

True

In [10]:

abs(gm(3.5) - mt.sqrt(np.pi) * 15 / 8) < 0.001

Out[10]:

True

In [11]:

abs(gm(4.5) - mt.sqrt(np.pi) * 105 / 16) < 0.001

Out[11]:

True

In [12]:

abs(gm(0.5) - mt.sqrt(np.pi)) < 0.001

Out[12]:

True

In [13]:

np.array([gm(1), gm(2), gm(3), gm(4), gm(5), gm(6), gm(7), gm(8)]).astype('int32')

Out[13]:

array([   1,    1,    2,    6,   24,  120,  720, 5040])

In [14]:

fact(np.arange(0,8)).astype('int32')

Out[14]:

array([   1,    1,    2,    6,   24,  120,  720, 5040])

In [15]:

### ### ###

Now let's explore the integrand function
$f_z(x) = x^{z-1}\cdot e^{-x}$
(which is used to define $\Gamma(z)$)
for several distinct fixed values of $z > 1$.

In this case we are exploring for values of $z$ such that $z \in [5, 6]$.

See the graphs/curves below. One curve corresponds to one particular value of $z$.
Thus for each value of $z$, we have that $\Gamma(z)$ is the area that is below the
graph of $f_z(x)$ but above the $X$ axis.

Note that if we explore $f_z(x)$ for values of $z < 1$ by plotting $f_z(x)$, the
graphs/curves then would look completely different from the graphs/curves shown here.

In [16]:

X_MIN = 0.001
X_MAX = 20

x = np.linspace(X_MIN, X_MAX, 200)
fig, ax = plt.subplots(figsize=(11,8))
ax.set(xlim=(X_MIN-0.5, X_MAX+5))

for z in np.arange(5.0, 6.01, 0.1):
    z = round(z, 3)
    y = np.power(x, z-1) * np.exp(-x)
    ax.plot(x, y, linewidth=2.0, label='$y = f_z(x)$ = x^(z-1).exp(-x) for z=' + str(z))

plt.xlabel('X')
plt.ylabel('Y', rotation=0)

plt.legend(loc='upper right', prop={'size': 12.5})
plt.show()