ベータ関数を組み合わせで表現する｜Python

ベータ関数は確率論・統計学で頻繁に登場する特殊関数だ。ベータ分布の正規化定数として現れ、ベイズ統計の共役事前分布でも重要な役割を果たす。Python の標準ライブラリにはベータ関数がないが、ガンマ関数と組み合わせを使って計算できる。

ベータ関数の定義

ベータ関数 $B (a, b)$ は次の積分で定義される。

$B (a, b) = \int_{0}^{1} t^{a - 1} (1 - t)^{b - 1} d t$

この積分は $a > 0$ , $b > 0$ で収束する。

ガンマ関数との関係

ベータ関数はガンマ関数で表現できる。

$B (a, b) = \frac{Γ ( a ) \cdot Γ ( b )}{Γ ( a + b )}$

この関係式を使って Python で実装する。

import math

def beta(a, b):
    """ベータ関数 B(a, b) を計算"""
    return math.gamma(a) * math.gamma(b) / math.gamma(a + b)

print(beta(2, 3))     # 0.08333333333333333（= 1/12）
print(beta(1, 1))     # 1.0
print(beta(0.5, 0.5)) # 3.141592653589793（= π）

$B (1/2, 1/2) = π$ という美しい結果が得られる。これはガンマ関数の性質 $Γ (1/2) = π$ から導かれる。

対数ベータ関数

大きな引数ではオーバーフローを避けるため、対数ベータ関数を使う。

import math

def log_beta(a, b):
    """対数ベータ関数 ln B(a, b) を計算"""
    return math.lgamma(a) + math.lgamma(b) - math.lgamma(a + b)

# 大きな値でもオーバーフローしない
print(log_beta(100, 200))  # -202.43614612888504
print(log_beta(1000, 2000))  # -2079.9413946292893

# 元の値が必要なら指数を取る
print(math.exp(log_beta(2, 3)))  # 0.08333333333333331

組み合わせとの関係

ベータ関数は二項係数（組み合わせ）とも関係がある。

$B (a, b) = \frac{1}{( a + b - 1 ) ( a - 1 a + b - 2 )}$

整数の場合、次の関係が成り立つ。

$B (m, n) = \frac{( m - 1 )! ( n - 1 )!}{( m + n - 1 )!} = \frac{1}{n ( m m + n - 1 )}$

import math

def beta_from_comb(m, n):
    """組み合わせを使ってベータ関数を計算（整数のみ）"""
    return 1 / (n * math.comb(m + n - 1, m))

# ガンマ関数版と比較
print(beta(3, 4))           # 0.016666666666666666
print(beta_from_comb(3, 4)) # 0.016666666666666666（一致）

ベータ関数の対称性

ベータ関数は引数を入れ替えても値が変わらない。

$B (a, b) = B (b, a)$

print(beta(2, 5))  # 0.03333333333333333
print(beta(5, 2))  # 0.03333333333333333（一致）

積分の定義から、 $t \to 1 - t$ と変数変換すれば対称性が示せる。

漸化式

ベータ関数は次の漸化式を満たす。

$B (a, b) = \frac{a - 1}{a + b - 1} B (a - 1, b) = \frac{b - 1}{a + b - 1} B (a, b - 1)$

a, b = 5, 3

# 漸化式の確認
lhs = beta(a, b)
rhs = (a - 1) / (a + b - 1) * beta(a - 1, b)
print(f'{lhs:.10f}')  # 0.0119047619
print(f'{rhs:.10f}')  # 0.0119047619（一致）

応用例：ベータ分布の確率密度関数

ベータ分布 $Beta (α, β)$ の確率密度関数は次のように定義される。

$f (x; α, β) = \frac{x ^{α - 1} ( 1 - x ) ^{β - 1}}{B ( α , β )}$

ベータ関数が正規化定数として現れる。

import math

def beta_pdf(x, alpha, beta_param):
    """ベータ分布の確率密度関数"""
    if x <= 0 or x >= 1:
        return 0
    numerator = x**(alpha - 1) * (1 - x)**(beta_param - 1)
    denominator = beta(alpha, beta_param)
    return numerator / denominator

# Beta(2, 5) の x = 0.3 での密度
print(beta_pdf(0.3, 2, 5))  # 2.016840000000001

ベイズ統計では、ベータ分布は二項分布の共役事前分布として使われる。

応用例：不完全ベータ関数と累積分布

不完全ベータ関数 $B_{x} (a, b)$ は積分の上限を $x$ に制限したもので、ベータ分布の累積分布関数に相当する。標準ライブラリにはないが、SciPy で計算できる。

from scipy.special import betainc

# 正則化不完全ベータ関数 I_x(a, b) = B_x(a, b) / B(a, b)
# これはベータ分布の累積分布関数と同じ
alpha, beta_param = 2, 5

# P(X ≤ 0.3) を計算
print(betainc(alpha, beta_param, 0.3))  # 0.58014800000000006

完全ベータ関数 B(a, b)

$0$ から $1$ までの積分。正規化定数として使う。

不完全ベータ関数

B_{x} (a, b)

$0$ から $x$ までの積分。累積分布関数の計算に使う。

SciPy のベータ関数

SciPy には専用の関数が用意されている。

from scipy.special import beta as sp_beta, betaln

# ベータ関数
print(sp_beta(2, 3))  # 0.08333333333333333

# 対数ベータ関数
print(betaln(100, 200))  # -202.43614612888506

大量の計算が必要な場合や、数値的な安定性が重要な場合は SciPy を使うとよい。