Beta, Incomplete Functions & Factorial Example

Demonstrates beta(), log_beta(), incomplete_gamma(), incomplete_gamma_upper(), incomplete_beta(), factorial(), and polygamma().

Example Code

"""Example: IMSL Special Functions — Beta, Incomplete Gamma/Beta, Factorial, Polygamma."""
from __future__ import annotations

import math
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import (
    beta, log_beta, incomplete_gamma, incomplete_gamma_upper,
    incomplete_beta, factorial, polygamma,
)

# ---------------------------------------------------------------------------
# Known-value checks
# ---------------------------------------------------------------------------
print(f"beta(2, 3)               = {beta(2, 3):.10f}  (expected: {1/12:.10f})")
print(f"log_beta(2, 3)           = {log_beta(2, 3):.10f}  (expected: {math.log(1/12):.10f})")
print(f"incomplete_gamma(1, 1)   = {incomplete_gamma(1, 1):.10f}  (expected: {1 - math.exp(-1):.10f}  = 1 - 1/e)")
print(f"incomplete_gamma_upper(1,1) = {incomplete_gamma_upper(1, 1):.10f}  (expected: {math.exp(-1):.10f})")
print(f"incomplete_beta(1, 1, 0.5)  = {incomplete_beta(1, 1, 0.5):.10f}  (expected: 0.5000000000)")
print(f"factorial(10)            = {factorial(10):.0f}  (expected: 3628800)")
print(f"factorial(0)             = {factorial(0):.0f}  (expected: 1)")
print(f"polygamma(0, 1)          = {polygamma(0, 1):.10f}  (expected: -0.5772156649  [digamma(1)])")
print(f"polygamma(1, 1)          = {polygamma(1, 1):.10f}  (expected: {math.pi**2/6:.10f}  [pi^2/6])")

# ---------------------------------------------------------------------------
# Plotting
# ---------------------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle("Beta, Incomplete Functions, Factorial & Polygamma", fontsize=14)

# --- Panel 1: Beta(a, b) vs a for several b values -------------------------
ax = axes[0, 0]
a_vals = np.linspace(0.2, 5.0, 400)
for b_val, colour in [(0.5, "tab:blue"), (1.0, "tab:orange"), (2.0, "tab:green")]:
    ax.plot(a_vals, beta(a_vals, b_val), color=colour, lw=1.8, label=f"b={b_val}")
ax.set_title("Beta(a, b) vs a")
ax.set_xlabel("a")
ax.set_ylabel("B(a, b)")
ax.set_ylim(0, 6)
ax.legend()
ax.grid(True, alpha=0.3)

# --- Panel 2: Regularised lower incomplete gamma P(a, x) vs x (Gamma CDF) --
ax = axes[0, 1]
x_vals = np.linspace(0, 10, 400)
for a_val, colour in [(1, "tab:blue"), (2, "tab:orange"), (3, "tab:green")]:
    ax.plot(x_vals, incomplete_gamma(a_val, x_vals), color=colour, lw=1.8,
            label=f"a={a_val}")
ax.set_title("Regularised Lower Incomplete Gamma P(a,x)\n(Gamma distribution CDF)")
ax.set_xlabel("x")
ax.set_ylabel("P(a, x)")
ax.legend()
ax.grid(True, alpha=0.3)

# --- Panel 3: Regularised incomplete beta I_x(a, b) vs x (Beta CDF) -------
ax = axes[1, 0]
x_vals_beta = np.linspace(0, 1, 400)
pairs = [(0.5, 0.5, "tab:blue"), (1, 1, "tab:orange"),
         (2, 5, "tab:green"), (5, 2, "tab:red")]
for a_val, b_val, colour in pairs:
    ax.plot(x_vals_beta, incomplete_beta(a_val, b_val, x_vals_beta),
            color=colour, lw=1.8, label=f"a={a_val}, b={b_val}")
ax.set_title("Regularised Incomplete Beta I_x(a,b)\n(Beta distribution CDF)")
ax.set_xlabel("x")
ax.set_ylabel("I_x(a, b)")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# --- Panel 4: Stirling approximation vs exact factorial (log scale) --------
ax = axes[1, 1]
ns = np.arange(1, 21)
exact = np.array([factorial(int(n)) for n in ns])
# Stirling: n! ≈ sqrt(2*pi*n) * (n/e)^n
stirling = np.sqrt(2 * np.pi * ns) * (ns / np.e) ** ns
ax.semilogy(ns, exact, "bo-", lw=1.5, ms=5, label="n! (exact)")
ax.semilogy(ns, stirling, "r--", lw=1.5, label="Stirling approx.")
ax.set_title("Stirling Approximation vs factorial(n)")
ax.set_xlabel("n")
ax.set_ylabel("n!  (log scale)")
ax.legend()
ax.grid(True, alpha=0.3, which="both")

plt.tight_layout()
plt.savefig("test_output/example_imsl_beta_incomplete.svg", bbox_inches="tight")
print("Saved: test_output/example_imsl_beta_incomplete.svg")

Plot Output

Beta(a,b) vs a, regularised lower incomplete Gamma P(a,x), regularised incomplete Beta I_x(a,b), and Stirling approximation compared to exact factorial(n).

Console Output

beta(2, 3)               = 0.0833333333  (expected: 0.0833333333)
log_beta(2, 3)           = -2.4849066498  (expected: -2.4849066498)
incomplete_gamma(1, 1)   = 0.6321205588  (expected: 0.6321205588  = 1 - 1/e)
incomplete_gamma_upper(1,1) = 0.3678794412  (expected: 0.3678794412)
incomplete_beta(1, 1, 0.5)  = 0.5000000000  (expected: 0.5000000000)
factorial(10)            = 3628800  (expected: 3628800)
factorial(0)             = 1  (expected: 1)
polygamma(0, 1)          = -0.5772156649  (expected: -0.5772156649  [digamma(1)])
polygamma(1, 1)          = 1.6449340668  (expected: 1.6449340668  [pi^2/6])
Saved: test_output/example_imsl_beta_incomplete.svg