ガンマ分布に関連した数値計算のあれこれ

ガンマ分布とは？

ガンマ分布 (Gamma distribution) とは指数分布を一般化した確率分布で

• 人の体重
• ウイルスの潜伏期間
• 電子部品の寿命
• トラフィックの待ち時間
• 所得

などがガンマ分布に従うと言われています。形状パラメータ $k > 0$、尺度パラメータ $\theta > 0$ を用いて、確率分布関数は

$$\mathit{Ga}(x | k, \theta) = \frac{1}{\Gamma(k) \theta^k} x^{k-1} e^{-x / \theta} \qquad (x > 0)$$

と表されます。ガンマ関数 $\Gamma(k)$ は階乗を連続量に拡張したもので、自然数について $k! = \Gamma(k+1)$ が成立します。

$$\Gamma(z) = \int^\infty_0 t^{z-1} e^{-t} \mathrm{d}t$$

ガンマ分布はパラメータ $\lambda = 1/\theta$ により

$$\mathit{Ga}(x | k, \lambda) = \frac{\lambda^k}{\Gamma(k)} x^{k-1} e^{-\lambda x}$$

と定義されることもあります。

ガンマ分布の予備知識

指数分布との関係性

ある期間に平均して $\lambda(= 1/\theta)$ 回起こる事象が、1回起きるまでの期間 $X_i$ が指数分布に従うとします。

$$X_i \sim \lambda e^{-\lambda x}$$

この確率変数が $k$ 個存在し、$X_1, X_2, \dots, X_k$ が独立な指数分布に従っていると考えます。このとき、これらの確率変数の和はガンマ分布に従います。

$$X_1 + X_2 + \cdots X_k \sim \mathit{Ga}(k, \lambda)$$

このため、ガンマ分布は「ある期間に平均して $\lambda$ 回起こる事象が $k$ 回起きるまでの期間」が従う分布と解釈できます。

※ $k$ が自然数の場合はアーラン分布であり、上記は厳密にはアーラン分布の説明です。ガンマ分布は $k$ を正の実数に拡張した分布なので、アーラン分布をさらに一般化した確率分布になっています。

平均と分散

• $E(X) = k \theta = k / \lambda$
• $V(X) = k \theta^2 = k / \lambda^2$

再生性

$X_1 \sim \mathit{Ga}(k_1, \theta), X_2 \sim \mathit{Ga}(k_2, \theta)$ であるとき、

$$X_1 + X_2 \sim \mathit{Ga}(k_1 + k_2, \theta)$$

となります。

パラメータ推定

モーメント法

モーメントを元にパラメータ推定を行う方法をモーメント法 (method of moment) といいます。モーメント法それ自体の説明は「モーメント法と最尤推定法の関連性についてメモ」がわかりやすいので、詳しく知りたい方はご覧ください。ここではガンマ分布のパラメータ推定の導出過程だけ書きます。

ガンマ分布の $n$ 次モーメントは

$$\begin{align*} E[X^n] & = \int^\infty_0 x^n \frac{1}{\Gamma(k) \theta^k} x^{k-1} e^{-x / \theta} \mathrm{d}x \\ & = \int^\infty_0 \frac{1}{\Gamma(k) \theta^k} x^{n+k-1} e^{-x / \theta} \mathrm{d}x \\ & = \int^\infty_0 \frac{k (k+1) \cdots (k+n)}{\Gamma(k+n)} \frac{\theta^n}{\theta^{n+k}} x^{n+k-1} e^{-x / \theta} \mathrm{d}x \\ & = k (k+1) \cdots (k+n) \theta^n \int^\infty_0 \frac{1}{\Gamma(k+n) \theta^{n+k}} x^{n+k-1} e^{-x / \theta} \mathrm{d}x \\ & = k (k+1) \cdots (k+n) \theta^n \cdot 1 \end{align*}$$

なので、標本平均を $\mu$、標本分散を $\sigma^2$ とすると、

• $E[X] = k \theta = \mu$
• $E[X^2] = k (k+1) \theta^2 = \sigma^2 + \mu^2$

となるはずです。これを $k$, $\theta$ について解くと

• $k = \mu^2 / \sigma^2$
• $\theta = \sigma^2 / \mu$

となり、パラメータを求められます。最尤推定の方が精度は良いですが、モーメント法でもそこそこの精度でパラメータを推定できます。

最尤推定

最尤推定 (maximum likelihood estimation) とは、尤度関数を最大化するようなパラメータを求める手法です。$n$ 個のデータ $D$ が与えられたとき、ガンマ分布の対数尤度関数は

$$\begin{align*} \ln L(D | k, \theta) & = \ln\left( \prod_i \mathit{Ga}(x_i | k, \theta) \right) = \ln\left( \prod_i \frac{1}{\Gamma(k) \theta^k} x^{k-1} e^{-x / \theta} \right) \\ & = (k-1) \sum_i \ln x_i - n \ln \Gamma(k) - n k\ln \theta - \frac{1}{\theta} \sum_i x_i \\ & = n (k-1) \overline{\ln x_i} - n \ln \Gamma(k) - n k\ln \theta - \frac{n \overline{x_i}}{\theta} \end{align*}$$

となり、$\ln L(D | k, \theta)$ を最大化するような $k, \theta$ を求めれば良いことになります。

局所近似最大化法

最適化の方法はいくつかありますが、ここでは Minka (2002) で提案された局所近似最大化法を紹介します。途中計算は「高速なガンマ分布の最尤推定法について - RPubs」を参照してもらうとして、以下の式を繰り返し計算すると、ガンマ分布のパラメータを求めることができます。

$$\begin{align*} k_0 & = \frac{0.5}{\ln \overline{x} - \overline{\ln x}} \\ \theta_0 & = \frac{\overline{x}}{k_0} \\ k_{i+1} & = \left( \frac{1}{k_i} + \frac{\overline{\ln x} - \psi(k_i) - \ln \overline{x} + \ln k_i}{k_i^2 (1 / k_i - \psi'(k_i))} \right)^{-1} \\ \theta_{i+1} & = \frac{\overline{x}}{k_{i+1}} \end{align*}$$

ただし、$\psi$ はディガンマ関数であり、ガンマ関数の対数を微分したものです。

$$\psi(z) = \frac{\mathrm{d} \ln \Gamma(z)}{\mathrm{d} z}$$

Python で書くとこんな感じになります。

import numpy as np
from scipy.special import gamma, loggamma, digamma, polygamma

def estimate_gamma_parameters(
        x: np.ndarray,
        tol: float = 1e-6,
        max_iter: int = 1000
    ) -> tuple[float, float, float]:
    """
    Estimate parameters of Gamma distribution.

    Returns
    -------
    k: float
        The shape parameter.
    theta: float
        The scale parameter.
    log_likelihood: float
        The log likelihood.
    """
    n = len(x)
    mean_x = np.mean(x)
    mean_log_x = np.mean(np.log(x))
    a = np.log(mean_x) - mean_log_x
    k = 0.5 / a
    theta = mean_x / k

    # ll = log likelihood
    ll0 = n * (k - 1) * mean_log_x - n * loggamma(k) - n * k * np.log(theta) - n * mean_x / theta
    for _ in range(max_iter):
        k = 1 / (1/k + (-a - digamma(k) + np.log(k)) / (k * k * (1/k - polygamma(1, k))))
        theta = mean_x / k
        ll1 = n * (k - 1) * mean_log_x - n * loggamma(k) - n * k * np.log(theta) - n * mean_x / theta
        if np.abs(ll1 - ll0) < tol:
            break
        ll0 = ll1

    return (k, theta, ll1)

累積分布関数（不完全ガンマ関数）

ガンマ関数の累積分布関数は正則化された第1種不完全ガンマ関数 (lower incomplete gamma) となります。

$$F(z) = \int^z_0 \mathit{Ga}(x | k, \theta) \mathrm{d} x = \frac{\gamma(k, z / \theta)}{\Gamma(k)}$$

第1種不完全ガンマ関数はガンマ関数の積分区間 $[0,\infty]$ の上限を変数化した関数です。

$$\gamma(a, x) = \int^x_0 t^{a-1} e^{-t} \mathrm{d} t$$

正則化第1種不完全ガンマ関数

不完全ガンマ関数をガンマ関数で割った関数を正則化不完全ガンマ関数 (regularized incomplete gamma) といいます。

正則化不完全ガンマ関数の数値計算

正則化不完全ガンマ関数を数値的に求める方法はいくつかありますが、ここではテイラー展開を使う方法を紹介します。

$$\frac{\gamma(a, x)}{\Gamma(a)} = x^a e^{-x} \sum^\infty_{k=0} \frac{x^k}{\Gamma(a+k)}$$

import numpy as np
from scipy.special import gamma

def rligamma(a: np.ndarray, x: np.ndarray, num_iter: int = 1000) -> np.ndarray:
    """
    Regularized lower incomplete gamma function by Taylor expansion.
    γ(a, x) = x^a * exp(-x) * sum(k=0..infty, x^k / Γ(a+k))
    """
    u = np.power(x, a) * np.exp(-x) / gamma(a + 1)
    s = u
    for k in range(1, num_iter):
        u = u * (x / (a + k))
        s += u
    return s

他の計算方法は 吉田 (2012) が詳しいです。

正則化不完全ガンマ関数の逆関数の数値計算

正則化不完全ガンマ関数が累積密度関数であるという前提に立つと、単調増加であることがわかります。単調増加ということは、二分法で解を求めることができます。

$\hat{y}, a$ を与えて $\hat{y} = \gamma(a, \hat{x}) / \Gamma(a)$ なる $\hat{x}$ を求める場合、以下のようなアルゴリズムで近似値を得られます。この方法はガンマ分布に限らず、累積分布関数を計算可能な任意の確率分布に対して適用できます。

1. 最初に適当な区間 $x_\mathrm{min}, x_\mathrm{max}$ を決める。
   • 経験的に $x_\mathrm{min} = 0, x_\mathrm{max} = 1.3 a + 40$ が良さそう。
2. 以下を繰り返す。
   1. 区間の中間点 $x_c = (x_\mathrm{min} + x_\mathrm{max}) / 2$ と中間点における $y_c = \gamma(a, x_c) / \Gamma(a)$ を計算する。
   2. $y_c$ と $\hat{y}$ が十分近い値なら $\hat{x}$ を出力して終了する。
   3. $y_c < \hat{y}$ ならば $\hat{x} \in [x_c, x_\mathrm{max}]$ なので、$x_\mathrm{min} \leftarrow x_c$ として次のループへ
   4. $y_c > \hat{y}$ ならば $\hat{x} \in [x_\mathrm{max}, x_c]$ なので、$x_\mathrm{max} \leftarrow x_c$ として次のループへ

初期値の $x_\mathrm{max} = 1.3 a + 40$ が天下り的に登場したので、説明します。$0.99999999 = \gamma(a, x) / \Gamma(a)$ となる $x$ を色々な $a$ に対して計算すると、以下のようなグラフになり、$1.3a + 40$ がその上限になっています。任意の $a$ について成り立つかどうか確認していないのですが、個人的に試した範囲ではうまく動作しています。もっと $y=1$ ギリギリの値でも逆関数を計算したい場合は、切片を増やすと良いと思います。

x=1.3a+40 は y=0.99999999 となる x の上限になっている

def inverse_rligamma(
        a: float,
        y: float,
        tol: float = 1e-6,
        max_iter: int = 1000,
        rligamma_kwargs = {}
    ) -> float:
    """
    The inverse of regularized lower incomplete gamma function.
    """

    xmin = 0.0
    xmax = 1.3 * a + 40

    for _ in range(max_iter):
        xc = (xmax + xmin) / 2
        yc = rligamma(a, xc, **rligamma_kwargs)
        if np.abs(y - yc) < tol:
            break
        elif yc < y:
            xmin = xc
        else:
            xmax = xc

    return xc

参考文献

• ガンマ分布 - Wikipedia
• 指数分布 - Wikipedia
• ガンマ分布の意味や性質と導出をわかりやすく解説 - 数理経済学的 SPECIAL PROJECT
• 高速なガンマ分布の最尤推定法について - RPubs
• 吉田 年雄. 漸化式を用いる不完全ガンマ関数 $\gamma(v,x)$ の数値計算法の誤差解析. 情報処理学会論文誌, 53(8), pp. 1954--1961. (2012)