Pythonプログラミング（ステップ７・補足・ガンマ関数）

このページでは、統計分布などにしばしば登場するガンマ関数について説明する。

積分で定義されるガンマ関数

積分を使って「新しい」関数を定義する例はたくさんある。 例えば、正規分布$\frac{1}{\sqrt{2\pi }}{e}^{-t^2/2}$ を元に、その累積分布関数 $$F(x)=\frac{1}{\sqrt{2\pi }}{\int }_{-\mathrm{\infty }}^x{e}^{-t^2/2}dt$$ を定義できるが、$F(x)$をこれ以上簡単な数式で表すことはできない。

これは、積分範囲を「可動」にして、関数を生成する例であるが、 被積分関数の中に「可動」な変数を入れて、関数を生成することもできる。 その代表的な例がガンマ関数で、以下のように定義される： $$\mathrm{\Gamma }(z)={\int }_0^{\mathrm{\infty }}{t}^{z-1}{e}^{-t}dt$$ 被積分関数の中身の$z$が「可動」箇所で、$t$で積分しても、その部分の自由度が残るため、積分全体としては$z$の関数となる。

階乗の一般化

では、ガンマ関数の値を計算してみよう。 $z=1$の場合は $$\mathrm{\Gamma }(1)={\int }_0^{\mathrm{\infty }}{e}^{-t}dt={[-e^{-t}]}_0^{\mathrm{\infty }}=1$$ $z=1/2$の場合は、$t=y^2$と変数変換すると $$\mathrm{\Gamma }(1/2)={\int }_0^{\mathrm{\infty }}{t}^{-1/2}{e}^{-t}dt=2{\int }_0^{\mathrm{\infty }}{e}^{-y^2}dy=\sqrt{\pi }$$ と求まる。 その他の一般の $z$ に対しては値を求めることが難しい。

しかしながら、$z$が自然数の場合と、半整数の場合は、具体的な値を得ることができる。 まず、部分積分($e^{-t}$を積分、$t^z$を微分）すると、 $$\begin{array}{rc}\mathrm{\Gamma }(z+1)=& {\int }_0^{\mathrm{\infty }}{t}^z{e}^{-t}dt={[-t^z{e}^{-t}]}_0^{\mathrm{\infty }}+{\int }_0^{\mathrm{\infty }}z{t}^{z-1}{e}^{-t}dt\\ & =z{\int }_0^{\mathrm{\infty }}{t}^{z-1}{e}^{-t}dt=z\mathrm{\Gamma }(z)\end{array}$$ という関係式（漸化式）が得られる。それを繰り返し適用すると、$n$が自然数の場合 $$\mathrm{\Gamma }(n+1)=n\mathrm{\Gamma }(n-1)=n(n-1)\mathrm{\Gamma }(n-2)=\cdots =n(n-1)(n-2)\cdots \mathrm{\Gamma }(1)=n!$$ となる。すなわち $\mathrm{\Gamma }(n+1)$ は $n!$ を与えることがわかった。 同様に、半整数に対しても $$\mathrm{\Gamma }\left(\frac{2n+1}{2}\right)=\frac{2n-1}{2}\frac{2n-3}{2}\cdots \mathrm{\Gamma }(1/2)$$ となるので、具体的な値を得ることができる。

ガンマ関数の近似：スターリングの式

統計学や物理学などでは、大きな$n$に対する階乗($n!$)の振る舞いを知る必要性が出てくる。 そこで、$n!$の良い見積を得るために、その積分表示である$\mathrm{\Gamma }(n+1)$を使ってみよう。 積分変数を$t=zs$と変換すると $$\begin{array}{r}\mathrm{\Gamma }(z+1)={\int }_0^{\mathrm{\infty }}{t}^z{e}^{-t}dt={\int }_0^{\mathrm{\infty }}(sz{)}^z{e}^{-sz}zds\\ =z^{z+1}{\int }_0^{\mathrm{\infty }}{s}^z{e}^{-sz}ds=z^{z+1}{\int }_0^{\mathrm{\infty }}{e}^{z(\log (s)-s)}ds\end{array}$$

ここで、$e$の「上」の $z(\log (s)-s)$ は上に凸な形をしており、微分を取ると、$s=1$で最大になることがわかる。 これを、$s=1$の周りでテイラー展開すると $$z(\log (s)-s)=-z-\frac{z}{2}(s-1{)}^2+\cdots$$ である。そこで、 $$\mathrm{\Gamma }(z+1)=z^{z+1}{e}^{-z}{\int }_0^{\mathrm{\infty }}{e}^{-\frac{z}{2}(s-1{)}^2+\cdots }ds$$ と書き直して、$z$が正の大きな値の場合の積分の振る舞いを考えてみる。

$z$が十分大きければ、被積分関数は$s=1$の周りで鋭いピークを持つので、その値は以下の式でよく近似できるであろう： $${\int }_0^{\mathrm{\infty }}{e}^{-\frac{z}{2}(s-1{)}^2+\cdots }ds\approx {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{-\frac{z}{2}(s-1{)}^2}ds=\sqrt{\frac{2\pi }{z}}$$ ピークの幅が鋭くなれば、「遠く」の様子は関係なくなるであろうから、積分範囲を$-\mathrm{\infty }$まで広げても結果は変わらないはずである。 また、展開の高次の項も無視できるとした（本来ならば誤差の評価を行うべきであるが、ここでは割愛）。

ここまでをまとめると、スターリングの式として知られる $$n!=\mathrm{\Gamma }(n+1)\approx n^{n+1}{e}^{-n}\sqrt{\frac{2\pi }{n}}=n^{n+1/2}{e}^{-n}\sqrt{2\pi }$$ が得られた。

Pythonでは、ガンマ関数は数学ライブラリとして用意されており、積分を行わずとも、import math しておけば、math.gamma( )で呼び出すことができる。以下に、ガンマ関数の値とスターリングの近似値のそれぞれ対数を計算し、比較してプロットするコードの例を示す。

# coding: utf-8

import matplotlib.pyplot as plt
import numpy as np
import math

xdata = np.arange(1,10,0.01)
n = xdata.shape[0]

log_gamma = np.empty(n)
log_stirling = np.empty(n)

for i in range(n):
    x = xdata[i]
    log_gamma[i] = math.log(math.gamma(x+1))
    log_stirling[i] = (x+1/2)*math.log(x) - x + 1/2*math.log(2*math.pi)

plt.title("GAMMA FUNCTION")
plt.plot(xdata,log_gamma, color=(1.0,0,0.0), linewidth=1.0)
plt.plot(xdata,log_stirling, color=(0.0,0,1.0), linewidth=1.0)
plt.xlim(0,10)
plt.xlabel('X')
plt.ylabel('LOG GAMMA(X)')
plt.grid(True)
plt.show()

ベータ関数

ガンマ関数の親戚で、統計学などでよく登場する二変数の関数にベータ関数がある。 ベータ関数はいろいろな定義（表現）が可能であるが、ひとまず、ガンマ関数との関係 $$B(x,y)=\frac{\mathrm{\Gamma }(x)\mathrm{\Gamma }(y)}{\mathrm{\Gamma }(x+y)}$$ を押さえておけば、Pythonのmathライブラリを使って、

math.gamma(x)*math.gamma(y)/math.gamma(x+y)

で計算が可能である。