バスがランダムにやってくる停留所で、3本目のバスが来るまでにどれくらい待つでしょうか。1本のバスの待ち時間が「指数分布」に従うことは知られていますが、では「3本ぶんの待ち時間の合計」はどんな分布になるのでしょう。これに答えるのがガンマ分布（Gamma distribution）です。
ガンマ分布は、正の実数 $x > 0$ を確率変数の範囲にもつ連続確率分布で、「いくつかの待ち時間の和」を表すのにぴったりの形をしています。実は指数分布もカイ二乗分布もガンマ分布の特別な場合で、これらを統一的に眺められる「大きな器」になっています。
ガンマ分布を理解すると、こんな場面で役に立ちます。
- 待ち時間・寿命のモデル化: 故障までの時間、窓口での処理時間、降雨量など「正の値をとり、右に裾を引く」量を柔軟に表せます。
- ベイズ統計: ポアソン分布のレート $\lambda$（単位時間あたりの発生回数）の共役事前分布として使われ、観測のたびにパラメータの足し算だけで信念を更新できます。
- 統計的検定の土台: カイ二乗分布はガンマ分布の特別な場合なので、適合度検定や分散の推定にもつながります。
本記事で扱う内容は次のとおりです。
- 「指数分布の和」という直感からガンマ分布を導入する
- 形状 $k$・尺度 $\theta$ による定義とガンマ関数 $\Gamma(\cdot)$ の役割を整理する
- 平均 $k\theta$・分散 $k\theta^2$ を1行ずつ導出する
- 指数分布の一般化・カイ二乗分布との関係・再生性を確かめる
- ポアソン率の共役事前分布としての使い方と最尤推定を導く
- Python でガンマ分布を可視化し、実データにフィットする
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 指数分布とは？平均や分散の導出や可視化で解説 — ガンマ分布の「1個ぶん」にあたる分布
- ポアソン分布の性質・期待値・分散について解説 — ガンマ分布が共役事前になる相手
- ベータ分布の完全ガイド — 同じく共役事前として使われる分布（混同しやすいので比較に）
ガンマ分布の直感 — 待ち時間をいくつ足すか
ガンマ分布の定義式は少し複雑に見えますが、その正体はとてもシンプルです。出発点は指数分布です。
「平均して $\theta$ 分に1回イベントが起きる」状況を考えましょう。停留所のバス、コールセンターへの電話、放射性原子の崩壊などがその例です。このとき、次の1回が起きるまでの待ち時間は平均 $\theta$ の指数分布に従います。
では、$k$ 回ぶんの待ち時間を合計したらどうなるでしょうか。3本目のバスが来るまでの時間は、「1本目までの待ち時間 + 2本目までの待ち時間 + 3本目までの待ち時間」の和です。1個ずつはバラつきますが、足し合わせると山のある分布になりそうです。この「平均 $\theta$ の指数待ち時間を $k$ 個足したもの」こそがガンマ分布です。
$$ X = E_1 + E_2 + \dots + E_k, \quad E_i \sim \mathrm{Exp}(\theta) \;\Longrightarrow\; X \sim \mathrm{Gamma}(k, \theta) $$
このイメージを押さえておくと、後で出てくる「平均 = $k\theta$」「分散 = $k\theta^2$」「再生性」がすべて自然に腑に落ちます。$k$ 個ぶんの待ち時間なら、平均も $k$ 倍になるのは当たり前ですよね。
上の図は、平均 $\theta=1$ の指数分布から取り出したサンプルを $k$ 個ずつ足し合わせ、その合計のヒストグラム（灰色）を理論曲線（赤）と重ねたものです。$k=1$ ではただの指数分布（右下がり）ですが、$k=2,4$ と足す個数を増やすほど、左端から立ち上がって山ができ、右に裾を引く形になっていきます。ヒストグラムと理論曲線がぴたりと一致しており、「指数分布の和 = ガンマ分布」が実際に成り立っていることが確認できます。
この直感を頭に入れたうえで、いよいよ正式な定義を見ていきましょう。
ガンマ分布の数学的定義
直感は「待ち時間の和」ですが、数式では確率密度関数（PDF）で定義します。まず定義式を示し、そのあとで各部品が何をしているかを丁寧に分解します。
ガンマ分布は、正の実数をとる2つのパラメータ、形状パラメータ（shape） $k > 0$ と 尺度パラメータ（scale） $\theta > 0$ を用いて、次のように定義されます。
$$ \begin{equation} \mathrm{Gamma}(x \mid k, \theta) = \frac{1}{\Gamma(k)\, \theta^{k}}\, x^{k-1} e^{-x/\theta}, \quad x > 0 \end{equation} $$
ここで $\Gamma(k)$ はガンマ関数で、次のセクションで詳しく説明します。各部品の役割は次のとおりです。
- $x^{k-1}$ は左端での立ち上がりを決めます。$k$ が大きいほど原点付近で密度が抑えられ、山が右に移動します。
- $e^{-x/\theta}$ は右の裾の減衰を決めます。$\theta$ が大きいほどゆっくり減り、分布が横に広がります。
- $\dfrac{1}{\Gamma(k)\theta^k}$ は、全体の積分を1にするための規格化定数です。
「shape（形状）」と「scale（尺度）」という名前のとおり、$k$ は山の形（とがり具合・対称性）を、$\theta$ は横方向の伸び縮みを担当します。
左の図は尺度 $\theta=2$ を固定して形状 $k$ を変えたもの、右の図は形状 $k=3$ を固定して尺度 $\theta$ を変えたものです。左図を見ると、$k=1$ では原点で最大の右下がり、$k$ を大きくするにつれて山ができて右へ移動し、なだらかになっていきます。右図では、$\theta$ を大きくすると山の位置が右へずれ、全体が横に広がって低くなります（面積は常に1なので、広がれば高さは下がります）。$k$ は「形」、$\theta$ は「スケール」という役割分担がはっきり見て取れます。
レート（rate）パラメータによる別表記
scipy や教科書によっては、尺度 $\theta$ の代わりにレートパラメータ $\beta = 1/\theta$ を使い、次のように書くこともあります。
$$ \mathrm{Gamma}(x \mid \alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)}\, x^{\alpha-1} e^{-\beta x} $$
ここで $\alpha = k$（形状）、$\beta = 1/\theta$（レート）です。ベイズ統計の文献ではこの $(\alpha, \beta)$ 表記が好まれます。意味は同じで、$\theta$ が「平均待ち時間」を表すのに対し、$\beta$ は「単位時間あたりのレート」を表します。本記事では直感に合う形状 $k$・尺度 $\theta$を主に使い、共役事前の節でレート表記にも触れます。
scipy の stats.gamma では形状を a、尺度を scale で指定します（scale=1/rate）。ライブラリを使うときは、自分がどちらの表記でパラメータを渡しているかを必ず確認しましょう。
定義が分かったところで、式の中に現れた $\Gamma(k)$ の正体を見ておきましょう。
ガンマ関数 — 階乗の連続版
定義式の分母にある $\Gamma(k)$ は、確率分布の規格化に欠かせないガンマ関数（gamma function）です。これは階乗 $n!$ を実数（さらには複素数）へ拡張したもので、次の積分で定義されます。
$$ \begin{equation} \Gamma(k) = \int_0^{\infty} t^{k-1} e^{-t}\, dt, \quad k > 0 \end{equation} $$
積分の中身 $t^{k-1}e^{-t}$ は、ガンマ分布の密度の「核」とそっくりです。実はガンマ関数は「ガンマ分布の核を全区間で積分した値」であり、だからこそ規格化定数に登場します。
ガンマ関数の最も重要な性質は、漸化式です。部分積分を1回行うと導けます。
$$ \Gamma(k+1) = k\,\Gamma(k) $$
これを繰り返し使うと、整数 $n$ に対して次が成り立ちます。
$$ \Gamma(n) = (n-1)! $$
つまりガンマ関数は階乗を「整数のすき間」まで滑らかにつないだものなのです。例えば $\Gamma(1)=0!=1$、$\Gamma(2)=1!=1$、$\Gamma(3)=2!=2$、$\Gamma(4)=3!=6$ となります。また半整数では $\Gamma(1/2)=\sqrt{\pi}$ という美しい値をとり、これがカイ二乗分布の係数に効いてきます。
青い曲線がガンマ関数 $\Gamma(x)$ で、赤い点は整数点での値 $\Gamma(n)=(n-1)!$ です。$x=1,2,3,4,5$ で $1, 1, 2, 6, 24$ を通っており、階乗の値をなめらかにつないでいることが分かります。$x$ が小さいところで $\Gamma(x)$ が大きくなるのは、原点付近で $\Gamma(x)\approx 1/x$ と発散するためです。ガンマ分布の密度が $1/\Gamma(k)$ で規格化されているのは、この曲線の値で割って面積を1に合わせているのだと理解できます。
ガンマ関数の道具がそろったので、ガンマ分布の期待値と分散を導出していきましょう。
ガンマ分布の期待値
ガンマ分布の期待値（平均）は、驚くほどシンプルな形になります。
$$ \begin{equation} E[X] = k\theta \end{equation} $$
「指数分布（平均 $\theta$）を $k$ 個足したもの」という直感を思い出せば、平均が $k$ 倍の $k\theta$ になるのは当然です。これを定義から確かめましょう。
ゴールは、密度に $x$ を掛けて積分し、$k\theta$ を得ることです。
$$ \begin{equation} \begin{split} E[X] &= \int_0^{\infty} x \cdot \frac{1}{\Gamma(k)\theta^k} x^{k-1} e^{-x/\theta}\, dx \\ &= \frac{1}{\Gamma(k)\theta^k} \int_0^{\infty} x^{k} e^{-x/\theta}\, dx \end{split} \end{equation} $$
ここで積分 $\int_0^{\infty} x^{k} e^{-x/\theta}\, dx$ を計算します。$t = x/\theta$ と置換すると $x = \theta t$、$dx = \theta\, dt$ なので、
$$ \int_0^{\infty} x^{k} e^{-x/\theta}\, dx = \int_0^{\infty} (\theta t)^{k} e^{-t}\, \theta\, dt = \theta^{k+1} \int_0^{\infty} t^{k} e^{-t}\, dt = \theta^{k+1}\, \Gamma(k+1) $$
最後はガンマ関数の定義式（指数が $k$ なので $\Gamma(k+1)$）を使いました。これを元の式に戻します。
$$ \begin{equation} \begin{split} E[X] &= \frac{1}{\Gamma(k)\theta^k} \cdot \theta^{k+1}\, \Gamma(k+1) \\ &= \theta \cdot \frac{\Gamma(k+1)}{\Gamma(k)} \end{split} \end{equation} $$
ここで漸化式 $\Gamma(k+1) = k\,\Gamma(k)$ を使うと $\Gamma(k+1)/\Gamma(k) = k$ なので、
$$ E[X] = \theta \cdot k = k\theta $$
無事に直感どおりの結果 $E[X]=k\theta$ が得られました。同じ要領で、分散も求められます。
ガンマ分布の分散
ガンマ分布の分散も簡潔な形になります。
$$ \begin{equation} V[X] = k\theta^2 \end{equation} $$
これも直感と整合します。独立な確率変数の和の分散は各分散の和なので、分散 $\theta^2$ の指数分布を $k$ 個足せば分散は $k\theta^2$ になります。定義から確かめましょう。分散の公式 $V[X] = E[X^2] – (E[X])^2$ を使うため、まず $E[X^2]$ を計算します。
$$ \begin{equation} \begin{split} E[X^2] &= \int_0^{\infty} x^2 \cdot \frac{1}{\Gamma(k)\theta^k} x^{k-1} e^{-x/\theta}\, dx \\ &= \frac{1}{\Gamma(k)\theta^k} \int_0^{\infty} x^{k+1} e^{-x/\theta}\, dx \end{split} \end{equation} $$
期待値のときと同じ置換 $t=x/\theta$ を使うと、今度は指数が $k+1$ なので $\Gamma(k+2)$ が現れます。
$$ \int_0^{\infty} x^{k+1} e^{-x/\theta}\, dx = \theta^{k+2}\, \Gamma(k+2) $$
これを戻して整理します。
$$ \begin{equation} \begin{split} E[X^2] &= \frac{1}{\Gamma(k)\theta^k} \cdot \theta^{k+2}\, \Gamma(k+2) \\ &= \theta^2 \cdot \frac{\Gamma(k+2)}{\Gamma(k)} \end{split} \end{equation} $$
漸化式を2回使うと $\Gamma(k+2) = (k+1)\Gamma(k+1) = (k+1)k\,\Gamma(k)$ なので、$\Gamma(k+2)/\Gamma(k) = k(k+1)$ です。よって、
$$ E[X^2] = \theta^2 \cdot k(k+1) = k(k+1)\theta^2 $$
最後に分散の公式へ代入します。$E[X]=k\theta$ より $(E[X])^2 = k^2\theta^2$ なので、
$$ \begin{equation} \begin{split} V[X] &= E[X^2] – (E[X])^2 \\ &= k(k+1)\theta^2 – k^2\theta^2 \\ &= (k^2 + k – k^2)\theta^2 \\ &= k\theta^2 \end{split} \end{equation} $$
こうして分散 $V[X]=k\theta^2$ が導けました。標準偏差は $\sqrt{k}\,\theta$ です。
上の図は、$(k,\theta)$ を変えた3つのガンマ分布について、平均 $k\theta$ を破線で、標準偏差 $\sqrt{k}\theta$ の幅を矢印で示したものです。破線が分布のおおよそ中心（やや左寄り、右裾があるため）に立っていること、$k$ が大きいほど（同じ平均でも）相対的なばらつき（変動係数 $1/\sqrt{k}$）が小さくなって分布が対称な山に近づくことが読み取れます。実際、$k\to\infty$ ではガンマ分布は正規分布に近づきます。
期待値と分散がそろったので、次はガンマ分布が他の有名な分布とどうつながっているかを見ていきます。
指数分布の一般化
冒頭の直感のとおり、ガンマ分布は指数分布を含む「一般化」です。定義式で $k=1$ とおいてみましょう。$\Gamma(1)=1$、$x^{1-1}=x^0=1$ なので、
$$ \mathrm{Gamma}(x \mid 1, \theta) = \frac{1}{\Gamma(1)\theta^1} x^{0} e^{-x/\theta} = \frac{1}{\theta} e^{-x/\theta} $$
これはまさに平均 $\theta$ の指数分布の密度です。つまり指数分布は「待ち時間1個ぶん」のガンマ分布 $\mathrm{Gamma}(1,\theta)$ にほかなりません。
赤い太線が $\mathrm{Gamma}(k=1,\theta=2)$、青い破線が指数分布 $\mathrm{Exp}(\theta=2)$ で、両者は完全に重なっています。薄い緑の曲線は $k$ を1より大きくしたガンマ分布で、$k=1$ の指数分布を「足し合わせて山を作った」ものだと分かります。指数分布は「メモリーレス（記憶がない）」という特別な性質をもちますが、$k>1$ のガンマ分布はその性質を失う代わりに、より柔軟に待ち時間の分布を表現できます。
指数分布が $k=1$ の特別な場合なら、もう一つ重要な特別な場合があります。カイ二乗分布です。
カイ二乗分布との関係
統計的検定で頻出するカイ二乗分布（chi-squared distribution）も、実はガンマ分布の特別な場合です。自由度 $n$ のカイ二乗分布は、次のガンマ分布に一致します。
$$ \chi^2_n = \mathrm{Gamma}\!\left(k = \frac{n}{2},\ \theta = 2\right) $$
確かめてみましょう。自由度 $n$ のカイ二乗分布の密度は次のとおりです。
$$ f_{\chi^2}(x) = \frac{1}{2^{n/2}\,\Gamma(n/2)}\, x^{n/2 – 1} e^{-x/2} $$
一方、ガンマ分布の定義式で $k=n/2$、$\theta=2$ を代入すると、
$$ \mathrm{Gamma}(x \mid n/2, 2) = \frac{1}{\Gamma(n/2)\, 2^{n/2}}\, x^{n/2 – 1} e^{-x/2} $$
両者は完全に一致します。カイ二乗分布は「標準正規分布に従う変数を2乗して $n$ 個足したもの」なので、これも「正の量の和」というガンマの仲間であることが効いています。
太線がカイ二乗分布、破線が $\mathrm{Gamma}(n/2, 2)$ で、各自由度 $n=2,4,6,9$ について両者がぴたりと重なっています。自由度が小さいほど右下がりに近く、自由度が大きいほど山ができて右に移動する様子は、ガンマ分布で $k=n/2$ を大きくしていくのと同じです。この関係を知っていると、カイ二乗分布の平均が $E[\chi^2_n] = (n/2)\cdot 2 = n$、分散が $(n/2)\cdot 2^2 = 2n$ であることが、ガンマの公式 $k\theta,\ k\theta^2$ から即座に分かります。
特別な場合を見たので、次はガンマ分布どうしを足したときの性質、再生性を確認します。
ガンマ分布の再生性
「待ち時間の和」という直感から予想できるとおり、尺度 $\theta$ が同じガンマ分布どうしの和は、再びガンマ分布になります。これを再生性（additivity / reproductive property）と呼びます。
$$ \begin{equation} X_1 \sim \mathrm{Gamma}(k_1, \theta),\quad X_2 \sim \mathrm{Gamma}(k_2, \theta) \;\Longrightarrow\; X_1 + X_2 \sim \mathrm{Gamma}(k_1 + k_2, \theta) \end{equation} $$
形状パラメータだけが足し算される点がポイントです。「$k_1$ 個ぶんの待ち時間」と「$k_2$ 個ぶんの待ち時間」を足せば「$k_1+k_2$ 個ぶんの待ち時間」になる、と思えば自然です。
この性質はモーメント母関数（MGF）を使うとすっきり証明できます。ガンマ分布のMGFは $t < 1/\theta$ に対して、
$$ M_X(t) = E[e^{tX}] = (1 – \theta t)^{-k} $$
です。$X_1, X_2$ が独立なら、和のMGFは積になります。
$$ M_{X_1 + X_2}(t) = M_{X_1}(t)\, M_{X_2}(t) = (1-\theta t)^{-k_1} (1-\theta t)^{-k_2} = (1-\theta t)^{-(k_1 + k_2)} $$
最後の式は $\mathrm{Gamma}(k_1+k_2, \theta)$ のMGFそのものです。MGFが一致すれば分布が一致するので、和は $\mathrm{Gamma}(k_1+k_2, \theta)$ に従います。
$X_1\sim\mathrm{Gamma}(2,1.5)$（青）と $X_2\sim\mathrm{Gamma}(3,1.5)$（緑）の和のヒストグラム（灰色）が、理論曲線 $\mathrm{Gamma}(5,1.5)$（赤）にぴたりと重なっています。形状が $2+3=5$ に足し算され、尺度 $\theta=1.5$ は変わらないことが実験的に確認できます。なお、尺度が異なるガンマ分布の和は一般にはガンマ分布になりません（待ち時間の「歩幅」が違うものは単純には足せない、とイメージするとよいでしょう）。
ここまでで分布としての性質が出そろいました。ここからは、ガンマ分布が活躍するベイズ統計の場面を見ていきましょう。
ポアソン率の共役事前分布
ガンマ分布がベイズ統計で重宝される最大の理由は、ポアソン分布のレート $\lambda$ の共役事前分布になることです。共役（conjugate）とは、「事前分布と事後分布が同じ分布族になる」性質のことで、これがあると計算が劇的に楽になります。
設定を整理しましょう。単位時間あたり平均 $\lambda$ 回起きるイベントを $n$ 回ぶん観測し、各区間での発生回数 $x_1, \dots, x_n$ が得られたとします。各 $x_i$ はポアソン分布に従います。
$$ p(x_i \mid \lambda) = \frac{\lambda^{x_i} e^{-\lambda}}{x_i!} $$
レート $\lambda$ の事前分布として、レート表記のガンマ分布 $\mathrm{Gamma}(\lambda \mid \alpha, \beta) \propto \lambda^{\alpha-1} e^{-\beta\lambda}$ を選びます。ベイズの定理より、事後分布は尤度と事前分布の積に比例します。
$$ \begin{equation} \begin{split} p(\lambda \mid x_{1:n}) &\propto \left[\prod_{i=1}^{n} \lambda^{x_i} e^{-\lambda}\right] \cdot \lambda^{\alpha-1} e^{-\beta\lambda} \\ &= \lambda^{\sum_i x_i} e^{-n\lambda} \cdot \lambda^{\alpha-1} e^{-\beta\lambda} \end{split} \end{equation} $$
ここで指数部分をまとめます。$\lambda$ の指数どうし、$e$ の肩どうしを足し合わせると、
$$ p(\lambda \mid x_{1:n}) \propto \lambda^{\,\alpha + \sum_i x_i – 1}\, e^{-(\beta + n)\lambda} $$
これは形状 $\alpha + \sum_i x_i$、レート $\beta + n$ のガンマ分布の核そのものです。したがって事後分布は、
$$ \begin{equation} p(\lambda \mid x_{1:n}) = \mathrm{Gamma}\!\left(\lambda \;\middle|\; \alpha + \sum_{i=1}^{n} x_i,\ \beta + n\right) \end{equation} $$
となり、事前と同じガンマ分布族に戻りました。これが共役性です。更新ルールは驚くほど単純で、「形状に観測総数 $\sum x_i$ を足し、レートに観測区間数 $n$ を足す」だけです。
上のフロー図は、事前分布 $\mathrm{Gamma}(\alpha,\beta)$ にポアソン観測を取り込むと、パラメータの足し算だけで事後分布 $\mathrm{Gamma}(\alpha+\sum x_i,\ \beta+n)$ が得られる様子を示しています。下のミニ図では、事前 $\mathrm{Gamma}(2,1)$（広く不確かな信念）が、データを取り込むことで事後 $\mathrm{Gamma}(20,6)$（鋭く絞られた信念）へと更新されています。データが増えるほど事後分布は真のレート付近に集中していきます。ちなみに、事後平均は $\dfrac{\alpha + \sum x_i}{\beta + n}$ となり、$n$ が大きくなると標本平均 $\frac{1}{n}\sum x_i$ に近づきます。
ベイズの枠組みを見たので、次はデータからパラメータを点推定する最尤推定を扱います。
ガンマ分布の最尤推定
観測データ $x_1, \dots, x_n$ から形状 $k$ と尺度 $\theta$ を推定するには、最尤推定（Maximum Likelihood Estimation, MLE）が基本です。対数尤度を最大化します。1サンプルの対数密度は、
$$ \log \mathrm{Gamma}(x \mid k, \theta) = (k-1)\log x – \frac{x}{\theta} – k\log\theta – \log\Gamma(k) $$
なので、$n$ 個ぶんの対数尤度は、
$$ \begin{equation} \ell(k, \theta) = (k-1)\sum_{i=1}^n \log x_i – \frac{1}{\theta}\sum_{i=1}^n x_i – n k \log\theta – n\log\Gamma(k) \end{equation} $$
まず $\theta$ について偏微分して0とおきます。
$$ \frac{\partial \ell}{\partial \theta} = \frac{1}{\theta^2}\sum_i x_i – \frac{nk}{\theta} = 0 $$
両辺に $\theta^2/n$ を掛けて整理すると、$\theta$ の推定値は形状 $k$ を使って次のように書けます。
$$ \hat\theta = \frac{1}{nk}\sum_{i=1}^n x_i = \frac{\bar x}{k} $$
ここで $\bar x$ は標本平均です。これを $\ell$ に代入して $k$ だけの関数（プロファイル対数尤度）にし、$k$ について偏微分して0とおくと、
$$ \log k – \psi(k) = \log \bar x – \overline{\log x} $$
という方程式が得られます。$\psi(k) = \frac{d}{dk}\log\Gamma(k)$ はディガンマ関数です。右辺は「平均の対数」と「対数の平均」の差で、データから直接計算できる定数です。この式は $k$ について閉じた解をもたないため、ニュートン法などで数値的に解きます。
左図はプロファイル対数尤度を $k$ の関数として描いたもので、1つのなめらかなピークをもち、その位置（赤破線）が最尤推定値 $\hat k$、真値（黒点線）とよく一致しています。右図は推定したガンマ分布をデータのヒストグラムに重ねたもので、フィット曲線（赤）が真の分布（黒点線）とほぼ重なり、データの形を的確に捉えていることが分かります。サンプル数が十分にあれば、最尤推定は真のパラメータをうまく回復してくれます。
理論がそろったので、ここからは Python で実際にガンマ分布を扱ってみましょう。
Pythonでの実装
新しい確率分布を学ぶときは、コードで触って動かしてみるのが一番です。scipy.stats.gamma を使い、可視化・統計量の確認・サンプリングを行います。
まず必要なライブラリを読み込み、形状 $k$ と尺度 $\theta$ を変えてガンマ分布を描いてみます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
x = np.linspace(0, 25, 500)
plt.figure(figsize=(10, 5))
# 尺度 theta=2 を固定して形状 k を変える
for k in [1, 2, 3, 5, 9]:
# scipy: a=形状, scale=尺度 (=1/レート)
y = stats.gamma(a=k, scale=2).pdf(x)
plt.plot(x, y, lw=2, label=f"k={k}, theta=2")
plt.xlabel("x")
plt.ylabel("密度 p(x)")
plt.legend()
plt.grid(alpha=0.3)
plt.show()
このコードを実行すると、$k=1$ では原点で最大の右下がり、$k$ を増やすと山ができて右へ移動するガンマ分布が描けます。a が形状、scale が尺度に対応している点に注意してください（レート $\beta$ を使いたいときは scale=1/beta と渡します）。本記事の冒頭の形状ギャラリーは、これと同じコードで生成しています。
次に、理論で導いた平均 $k\theta$・分散 $k\theta^2$ が、scipy の計算やモンテカルロと一致するか確認します。
k, theta = 3.0, 2.0
dist = stats.gamma(a=k, scale=theta)
# 理論値
print(f"理論 平均 = k*theta = {k*theta:.4f}")
print(f"理論 分散 = k*theta^2 = {k*theta**2:.4f}")
# scipyの計算
print(f"scipy 平均 = {dist.mean():.4f}")
print(f"scipy 分散 = {dist.var():.4f}")
# モンテカルロ（サンプリングして標本統計量）
rng = np.random.default_rng(0)
samples = rng.gamma(shape=k, scale=theta, size=200000)
print(f"標本 平均 = {samples.mean():.4f}")
print(f"標本 分散 = {samples.var():.4f}")
実行すると、理論値・scipy・標本値がいずれも平均 $\approx 6.0$、分散 $\approx 12.0$ で一致します（$k\theta = 6$、$k\theta^2 = 12$）。導出した公式が正しいことを数値で裏づけられました。標本統計量はサンプル数を増やすほど理論値に近づきます。
指数分布の和がガンマになることの確認
冒頭の直感「指数分布を $k$ 個足すとガンマ分布」を、実験で確かめてみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(1)
theta = 1.0
k = 4 # 足し合わせる指数分布の個数
N = 60000
# 平均theta の指数分布を k個サンプリングして合計する
expo = rng.exponential(scale=theta, size=(N, k))
sums = expo.sum(axis=1)
x = np.linspace(0, 12, 400)
plt.figure(figsize=(9, 5))
plt.hist(sums, bins=80, density=True, alpha=0.5, label=f"指数{k}個の和（実験）")
plt.plot(x, stats.gamma(a=k, scale=theta).pdf(x), "r", lw=2.5,
label=f"Gamma(k={k}, theta={theta})")
plt.xlabel("待ち時間の合計"); plt.ylabel("密度")
plt.legend(); plt.grid(alpha=0.3); plt.show()
ヒストグラム（実験）と理論曲線 $\mathrm{Gamma}(4, 1)$ がぴたりと重なります。これで「ガンマ分布 = 指数待ち時間の和」が単なる比喩ではなく、数学的に成り立つ事実だと確認できました。k を変えると山の位置や形が変わることもぜひ試してみてください。
実データへのフィット
最後に、実データ（ここでは窓口の処理時間を模した合成データ）にガンマ分布を最尤推定でフィットします。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
rng = np.random.default_rng(2024)
# 真の分布: Gamma(k=2.5, theta=4.0) からの観測（待ち時間データの想定）
data = rng.gamma(shape=2.5, scale=4.0, size=600)
# scipyで最尤推定（loc=0に固定して形状kと尺度thetaを推定）
k_hat, loc, theta_hat = stats.gamma.fit(data, floc=0)
print(f"推定 k = {k_hat:.3f}")
print(f"推定 theta = {theta_hat:.3f}")
print(f"推定 平均 = {k_hat*theta_hat:.3f} (標本平均 {data.mean():.3f})")
xx = np.linspace(0, data.max(), 400)
plt.figure(figsize=(9, 5))
plt.hist(data, bins=40, density=True, alpha=0.6, label="待ち時間データ")
plt.plot(xx, stats.gamma(a=k_hat, scale=theta_hat).pdf(xx), "r", lw=2.5,
label=f"ガンマフィット k={k_hat:.2f}, theta={theta_hat:.2f}")
plt.xlabel("待ち時間 [分]"); plt.ylabel("密度")
plt.legend(); plt.grid(alpha=0.3); plt.show()
推定された $\hat k \approx 2.5$、$\hat\theta \approx 4.0$ は真の値をよく回復し、フィット曲線がヒストグラムに沿って立ち上がりから裾まで捉えています。指数分布（$k=1$ に固定）では原点で最大の右下がりしか表せませんが、ガンマ分布は $k>1$ によって「最初は少なく、ピークを過ぎてから減る」という現実的な待ち時間の形を表現できるのが強みです。
左図がフィット結果で、ガンマフィット（赤）はデータの立ち上がりを捉えていますが、指数フィット（緑破線、$k=1$ 固定）は原点で最大になってしまい形が合いません。右図はガンマ分布のQ-Qプロットで、点がほぼ直線上に並んでいることから、データがガンマ分布によく適合していると判断できます。Q-Qプロットは「分布の仮定が妥当か」を視覚的にチェックする定番の方法です。
ガンマ分布の応用
ここまでの性質を踏まえて、ガンマ分布が実務でどう使われるかを整理しておきましょう。
- 待ち時間・寿命のモデル化: 機械が故障するまでの時間、サービス窓口での処理時間、地震や事故の発生間隔など、「正の値をとり右に裾を引く」量に広く使われます。$k$ を調整することで、指数分布では表せない「最初は起きにくく、ある程度たってから起きやすくなる」挙動も表現できます。信頼性工学のワイブル分布と並ぶ基本的な寿命分布です。
- 降雨量・水文学: 日降水量や河川流量はガンマ分布でよく近似されます。気象データの解析や、標準化降水指数（SPI）による干ばつ判定にもガンマ分布が使われています。
- ベイズ統計: ポアソン率の共役事前分布として、コールセンターの着信率、ウェブサイトのアクセス率、保険のクレーム発生率などの推定に使われます。また、正規分布の精度（分散の逆数）の共役事前分布もガンマ分布で、ガウス分布のベイズ推定に欠かせません（分散そのものには逆ガンマ分布を使います）。
- 金融・保険: 損害額や保険金請求額の分布のモデル化にも用いられます。
これらの応用に共通するのは、「正の量で、平均のまわりに非対称（右に裾）なばらつきをもつ」という特徴です。ガンマ分布の2つのパラメータ $k, \theta$ が、この多様な現象を柔軟に表現してくれます。
まとめ
本記事では、ガンマ分布について直感・定義・導出・応用まで解説しました。
- ガンマ分布は「平均 $\theta$ の指数待ち時間を $k$ 個足したもの」と理解できる、正の実数上の連続分布
- 形状 $k$ と尺度 $\theta$（またはレート $\beta=1/\theta$）の2パラメータで、規格化にはガンマ関数 $\Gamma(k)$ が現れる
- 平均は $k\theta$、分散は $k\theta^2$、標準偏差は $\sqrt{k}\theta$
- $k=1$ で指数分布、$k=n/2,\ \theta=2$ でカイ二乗分布に一致する一般化
- 尺度が同じガンマの和は形状が足し算されるガンマになる（再生性）
- ポアソン率 $\lambda$ の共役事前分布で、更新は「形状に総数、レートに区間数を足す」だけ
ガンマ分布は、指数分布・ポアソン分布・カイ二乗分布・逆ガンマ分布・ベータ分布など、多くの分布を結ぶハブのような存在です。ここを押さえておくと、確率分布の地図がぐっと見通しよくなります。
次のステップとして、以下の記事も参考にしてください。
- 指数分布とは？平均や分散の導出や可視化で解説 — ガンマ分布の $k=1$ の場合
- χ²分布の導出と性質 — 正規分布の二乗和が従う確率分布 — ガンマ分布の $k=n/2,\theta=2$ の場合
- 逆ガンマ分布完全ガイド — 共役事前分布としてのベイズ推定への応用 — 正規分布の分散の共役事前
- ベータ分布の完全ガイド — 二項分布の共役事前（混同しやすい分布）