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

Gamma function

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

for $x,t\in \mathbb R$ and $\Gamma: \mathbb R \to \mathbb R$. However, one can also extend the definition to complex numbers by analytic continuation and obtain the same definition but with $x, \Gamma(x) \in \mathbb C$. In Scipy, one uses the complex definition.

x = np.linspace(-2,6,1000)
G = gamma(x)

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

(-100.0, 100.0)

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

/var/folders/2m/tgbjydr93_9_19gzzmp_v_3w0000gn/T/ipykernel_18975/2765431682.py:3: RuntimeWarning: invalid value encountered in log
  plt.plot(x,np.log(gamma(x)))

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

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

/var/folders/2m/tgbjydr93_9_19gzzmp_v_3w0000gn/T/ipykernel_18975/3274935164.py:3: RuntimeWarning: invalid value encountered in log10
  plt.plot(x,np.log10(gamma(x)))

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

Comparison with the factorial:

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

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

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

24 24

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

x_posinteger = np.linspace(0,6,6,dtype=int)
n = x_posinteger + 1
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()
plt.ylim(-100,200)

(-100.0, 200.0)

Recursion of the Gamma function

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)

The gamma function and 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

Let's explore numerically

x = np.linspace(-10,10,100)
I0_integrand = np.exp(-x**2)
plt.plot(x,I0_integrand)

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

dx = np.diff(x)[0]
np.sum(I0_integrand*dx), np.sqrt(np.pi)

(np.float64(1.7724538509055194), np.float64(1.7724538509055159))

Beta function

from scipy.special import beta

beta(2,3)

np.float64(0.08333333333333333)

gamma(2)*gamma(3)/gamma(2+3)

np.float64(0.08333333333333333)

beta(3,2)

np.float64(0.08333333333333333)

print(beta(1,5), 1/5)

0.2 0.2

Complex gamma function

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)

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

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