import numpy as np
from scipy.special import gamma
import matplotlib.pyplot as plt

Gamma function

A generalization of the factorial

Let's remember the factorial. The factorial of an integer $n$ is defined recursively as

$$ n! = n(n-1)! $$

This recursion continues without end (we keep subtracting 1) until we reach 0!. We define $0! \equiv 1$. That is,

$$ n! = n\cdot (n-1) \cdot (n-2) \cdots 1 $$

In Python, we can program this with a recursive function:

def factorial(n):
    if n > 0:
        return n*factorial(n-1)
    if n == 0:
        return 1
    else:
        raise ValueError

Let's test it

print(factorial(4),4*3*2*1)

24 24

Now let's plot the factorials of integer numbers:

n = np.linspace(1,6,6,dtype=int)
factorials = np.array([ factorial(n_i) for n_i in n ])

plt.scatter(n,factorials)

<matplotlib.collections.PathCollection at 0x112ce4b00>

Now let's imagine there were a curve that is able to connect the dots. One such curve is called the Gamma function and it is defined recursively in a similar way than the factorial:

$$ \Gamma(x+1) = x \Gamma (x) $$

$\Gamma(x) : \mathbb R \to \mathbb R$.

We can relate this definition to the factorial for integer $n$ as

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

Please don't forget the -1 when converting from factorial to gamma and vice versa!

Let's take a look at what this function looks like for $x > 0$ and the comparison to the factorial function, using Scipy.

x = np.linspace(0,6,1000)
G = gamma(x+1)

x_posinteger = np.linspace(0,6,6,dtype=int)
n = x_posinteger
factorials = np.array([ factorial(n_i) for n_i in n ])

plt.plot(x,G,label="Gamma(x+1)")
plt.scatter(n,factorials,label="n!")
plt.legend()

<matplotlib.legend.Legend at 0x112f0dfd0>

How to define the Gamma function?

The Gamma function is defined as

$$ \Gamma(x) = \int_0^\infty t^{x-1} e^{-t} dt $$

$\Gamma(x): \mathbb R \to \mathbb R$, $t\in \mathbb R$.

Important analytical exercise: show the properties of a) recursion b) $\Gamma(1) = 1$ from this definition.

(a) Recursivity: we have $\Gamma(x+1) = \int_0^\infty t^{x} e^{-t} \, dt $. Integrating by parts, we have $$ \Gamma(x+1) = \left[-t^x e^{-t}\right]_0^\infty + \int_0^\infty x \, t^{x-1} e^{-t}\, dt = x \Gamma(x)$$ after seeing that the evaluation at the limits goes to zero and using the definition of Gamma.

(b) $\Gamma(1) = \int_0^\infty t^{1-1} e^{-t}\,dt = \int_0^\infty e^{-t} \, dt = 1.$

Let's see how the Gamma function then looks like for $x \in \mathbb R$.

x = np.linspace(-3,6,5000)
G = gamma(x)

plt.xticks(range(-3,7))

plt.plot(x,G)
plt.ylim(-100,100)

(-100.0, 100.0)

We see that Gamma of negative integers [e.g. $\Gamma(-2)$] is not defined, but Gamma of negative non-integers [e.g. $\Gamma(-\pi/2)$] is defined.

Recursivity is part of the definition of $\Gamma(x)$. In this plot, we test that the function is indeed recursive.

x = np.linspace(-2,6,1000)
G = gamma(x+1)
G1 = x*gamma(x)
plt.plot(x,G,ls="-",label="Gamma(x+1)")
plt.plot(x,G1,ls="--",label="x*Gamma(x)")
plt.legend()
plt.ylim(-100,100)

(-100.0, 100.0)

What happens when the argument is a large number?

In the plot for positive arguments of the Gamma function, one can suspect that there is some sort of nice exponential behavior. First, let's see again what happens to larger and larger factorials:

n = np.linspace(1,8,8,dtype=int)
factorials = np.array([ factorial(n_i) for n_i in n ])

plt.scatter(n,factorials)

<matplotlib.collections.PathCollection at 0x113066d20>

The order of magnitude of the factorial increases with each new number in the horizontal axis. Let's try to visualize this by taking the logarithm base 10 of the Gamma function:

# Logarithm base 10 of the Gamma function
x = np.linspace(0,20,1000)
plt.plot(x,np.log10(gamma(x)))
# easy way to visualize n! of a very large number

[<matplotlib.lines.Line2D at 0x1131f17c0>]

Taking the logarithm has made the input of large numbers manageable. The natural logarithm also works:

# Natural logarithm of the Gamma function
x = np.linspace(0,20,1000)
plt.plot(x,np.log(gamma(x)))

[<matplotlib.lines.Line2D at 0x1132f5a30>]

Important analytical exercise: because of this very nice behavior for large numbers, one can approximate the factorial of large numbers (or equivalently, the Gamma function of large numbers) as

$$ \ln n! \approx n \ln n - n $$

This is called the Stirling approximation. We can visualize it in the following plot.

n = np.linspace(2,50,1000)
plt.plot(n+1,np.log(gamma(n+1)), label=r"$\ln(n!)$")
plt.plot(n,n*np.log(n)-n, label=r"Stirling's formula")
plt.legend()

<matplotlib.legend.Legend at 0x1131f20f0>

The Gamma function and the Gaussian distribution

Important analytical exercise: In Physics, we frequently encounter the Gaussian distribution, $\sim e^{-x^2}$. Evaluate $\Gamma(1/2)$ by using the Gaussian distribution.

Consider the integral $I_0 = \int_{-\infty}^{\infty} e^{-a x^2} dx$. We can integrate by considering first the integral $I_0^2 = \int_{-\infty}^\infty e^{a(x^2 + y^2)} dx dy$ and changing to polar coordinates, $I_0^2 = \int_0^{\infty} \int_0^{2\pi} e^{-ar^2} rdrd\phi = \pi/a $ $\implies I_0 = \sqrt{\pi/a}$.

For $a = 1$, we have $I_0 = \sqrt{\pi}$. Now think of the integral

$$ \Gamma(1/2) = \int_0^{\infty} \frac{1}{\sqrt{t}}e^{-t} dt $$

With the change of variable $t = y^2$, we find $ \Gamma(1/2) = 2\int_0^\infty e^{-y^2} dy = \sqrt{\pi} $

# Verification
print(gamma(1/2),np.sqrt(np.pi))

1.7724538509055159 1.7724538509055159

By the recurrence relation, we can also find

print(gamma(3/2), (1/2)*gamma(1/2), (1/2)*np.sqrt(np.pi))

0.8862269254527579 0.8862269254527579 0.8862269254527579

What happens for complex numbers?

We can define $\Gamma(z): \mathbb C \to \mathbb C$, $t\in \mathbb R$, as

$$ \Gamma(z) = \int_0^\infty t^{z-1} e^{-t} dt $$

in a similar way than the real $\Gamma(x)$, through the integral definition.

x = np.linspace(-10,10,1000)
y = np.linspace(-10,10,1000)
x, y = np.meshgrid(x,y,indexing="ij")
z = x + 1j*y
gamma_eval = gamma(z)

Let's plot the real part of $\Gamma(z)$:

plt.pcolormesh(np.real(z),np.imag(z),np.real(gamma_eval),vmin=-10,vmax=10,cmap="seismic_r")
plt.colorbar()
plt.grid()

Mental check: in the horizontal axis, the function grows strongly just as we have plotted before (following the factorial function). For negative numbers, we see that the function becomes undefined for negative integers.

Now, let's plot the imaginary part:

plt.pcolormesh(np.real(z),np.imag(z),np.imag(gamma_eval),vmin=-10,vmax=10,cmap="seismic_r")
plt.colorbar()
plt.grid()

Mental check: along the horizontal axis, the imaginary part is everywhere zero, meaning that $\Gamma(x)$ is indeed real for real $x$.