確率分布には様々の特殊関数が用いられています．なかでもガンマ関数・ベータ関数などは必須項目であり，さらに非心分布ではベッセル関数・超幾何関数などの理解も必要になります．しかし，これらの特殊関数はプロの統計研究者にとっても馴染みの薄いものであり，ましてや一般のひとにとって（無論，私にとっても）とっつきにくいものです．

　ここでは，特殊関数に加え，ラプラス変換やメリン変換など数理統計学における方法的基礎を，いろいろな話題を交えながら，わかりやすく解説しきたいと思います．なお，副読本としては「特殊関数とその応用」（奥井著：森北出版）をお勧めします．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

(1)ガウス積分

　正規分布は一般的な誤差の分布関数で，その確率密度関数，累積分布関数それぞれ

f(x)=1/√2πexp(-x^2/2)=φ(x)

F(x)=∫(-∞,x)f(t)dt=Φ(x)

と表されます．ここでは，正規分布の累積分布関数Φ(x)に関連して，

I=∫(0,∞)exp(-x^2)dx

の値を計算してみます．

　ケルビン卿の銘言に「数学者とは

　　∫(-∞,∞)exp(-x2)dx=√π

を１＋１＝２のように自明だと思っている人である」とあります．

　われわれは数学者ではないですが，以下のように極座標を用いることによって，簡単に数学者になることができます．

I^2=∫(0,∞)exp(-x2)dx∫(0,∞)exp(-y2)dy（２重積分）

=∫(0,π/2)int(0,∞)exp(-r2)rdrdθ（極座標変換）

より，結局，I=√π/2となります．

　以前より，どうして正規分布に円周率πが現れるか疑問視しておられた方も多いと思いますが，極座標に変換することによって，πが自然に入り込んできます．また，ここでは２重積分を用いてガウス積分を解きましたが，複素積分を用いると，もっと直接的に角度と関係していることが理解されます．ともあれ，πは幾何のみならず，統計にも使われることになります．

　また，上式において，x^2=tとおくと，ガウス積分とガンマ関数との面白い関係

　　√π=2I=2∫(0,∞)exp(-x2)dx=∫(0,∞)exp(-t)/√tdt=Γ(1/2)

も得られます．

　なお，逆関数Φ-1(x)については

∫(0,1)Φ-1(x)dx=0

∫(0,1)[Φ-1(x)]^2dx=1

∫(0,1)xΦ-1(x)dx=1/(2√π)

が成り立ちます．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

(2)ガンマ関数

Γ(x)=∫(0,∞)t^(x-1)e^-tdt x>0

この無限積分をｘの関数とみてガンマ関数Γ(x)といいます．

Γ(1)=∫(0,∞)e^-tdt=1

Γ(1/2)＝∫(0,∞)t^(-1/2)e^-tdt

　ここで，t=u^2とおくと∫(0,∞)e^-u2/2du=√π/2（ガウス積分）より

　　Γ(1/2)＝√π

が得られます．

　オイラーの第２種積分とも呼ばれるガンマ関数Γ（ｘ）には，Γ（ｘ＋１）＝ｘΓ（ｘ）の関係があり，次のような漸化式が成り立ちます．

　　Γ(x+1)=xΓ(x)=x(x-1)Γ(x-1)=････

したがって，ｘが正の整数ｎのときには，

　　Γ（ｎ＋１）＝ｎ！

が成り立ち，ガンマ関数は階乗の一般形となっていることがわかります．階乗の解析的補間をしている関数がガンマ関数なのです．

（ガンマ関数と超球との関係）

ガウス積分をｎ次元に拡張し，

I=∫(-∞,∞)exp(-x1^2+x2^2+･･･+xn^2)dx1dx2･･･dxn

を考えると∫(-∞,∞)exp(-x2)dx=√πのｎ重積分より，直ちに

　　I=π^(n/2)

を得ることができます．

　n次元ガウス積分を別の方法，直交座標でなく，極座標で求めてみましょう．

球に相当するｎ次元の図形を超球と呼びます．ｎ次元単位超球{x12+x22+･･･+xn2≦1}の体積をＶnとすると，Ｖ1=2（直径）,Ｖ2=π（面積）,Ｖ3=4π/3（体積）はご存知でしょう．また，単位超球の表面積Ｓn-1はｎＶn，半径ｒのｎ次元球の体積はｖnｒ^n，表面積はｎＶnｒ^n-1となります．

　ガウス積分の被積分関数を原点を中心とする半径ｒの球面上で積分し，次にｒ＝０からｒ＝∞まで積分すると，半径ｒの球面上で被積分関数は一定値exp(-r^2)をとり，表面積はｎＶnｒ^n-1ですから，

I=∫(0,∞)exp(-r2)ｎＶnｒ^n-1dr=ｎＶn∫(0,∞)r^(n-1)exp(-r2)dr

z=r2と変数変換するとdz=2rdrより

I=ｎＶn/2∫(0,∞)z^(n/2-1)exp(-z)dz

=Ｖnn/2Γ(n/2) n/2Γ(n/2)=Γ(n/2+1)

=ＶnΓ(n/2+1)

　したがって，

　　Ｖn=π^(n/2)/Γ(n/2+1)　　　ｒn=

を得ることができます．この結果は，形式的にＶn=π^(n/2)/(n/2)！と書くことができます． Γ(m+1)=m!

これより，半径ｒのｎ次元超球の超体積はＶnｒn=(πr2)^(n/2)/Γ(n/2+1)となります．

　ｎが整数のとき，実際にＶnの値を計算してみると，超球の体積はｎ＝５のとき最大８π2／１５＝５．２６３７・・・となり，以後は減少します．

ｎ　　１次元 ２次元 ３次元　４次元　５次元　　６次元

Ｖn　　２ 　 3.14 4.19 　 4.93 5.263 5.167

ｎ　　７次元 ８次元 ９次元　１０次元

　　 4.72 　 4.06　　 3.30 　　2.55

（次元を整数に限らなければ５．２５６次元で最大となり，そのときの体積は5.277･･･である．）

　Ｖn-1がわかればＶnは漸化式：

Ｖn/Ｖn-1=Γ（1／２）Γ｛（ｎ+１）／２｝／Γ（ｎ／２+1）=B(1/2,(n+1)/2)

によって求めることができますが，この計算は面倒ですから，Ｖn-2との漸化式

　　Ｖn/Ｖn-2=2π/n

を用いると任意のｎに対して

ｎが奇数であればＶn=2(2π)^((n-1)/2)/n!!

ｎが偶数であればＶn=(2π)^(n/2)/n!!

とも書けることも理解されます．

　ｎ→∞のとき

　　Ｖn/Ｖn-2=2π/n→０，Ｓn-1/Ｓn-3=nＶn/(n-2)Ｖn-2=2π/(n-2)→０

ですから，不思議なことに，単位球面の体積や表面積はｎ→∞のとき０に収束するのです．

　また，このことから，ｎ次元単位超立方体[-1,1]^nにおいて，単位超球が占める比率は，ｎ＝２であればπ／４(79%)であるが，ｎ＝５のときは16%に下落し，ｎ＝１０となると0.25%になることも理解されます．高次元において，超立方体内に一様分布する標本を考えるとき，低次元の場合とは対照的に，大部分のデータは超球外に位置することになります．

　なお，比Ｖn-1/Ｖn-2=B(1/2,n/2)は自由度ｎのｔ分布の定数であり，実際，フィッシャーはｎ個の観測値の標本平均と母平均の差（距離）を標本標準偏差で割った統計量ｔの分布をｎ次元ユークリッド空間を使って導きだしています．

【補】ウォリスの公式

1/2B(1/2,(n+1)/2)=∫(0-π/2)(sinθ)^ndθ

この値をＳnとおくと，部分積分により漸化式

　　Ｓn=(n-1)/nＳn-2

が得られますから，

n=2k（偶数）なら1･3･･･(2k-1)/2･4･･･(2k)*π/2

n=2k+1（奇数）なら2･4･･･(2k)/1･3･･･(2k+1)

　これより，ｌｉｍ1･3･･･(2k-1)/2･4･･･(2k)*√(k)=1/√(π)

変形するとウォリスの公式

　　(2n)!/(2^nn!)^2√(n)=1/√(π)

が得られる．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

(3)ポリガンマ関数

　ガンマ関数の対数微分であるジガンマ関数φ（ｘ）は

φ(x)=d/dx{logΓ(x)}=Γ'(x)/Γ(x)

で定義されます．また，その逐次導関数φ’（ｘ），φ”（ｘ），・・・，φ^(k)（ｘ），すなわち，トリガンマ関数，テトラガンマ関数，ペンタガンマ関数などを総称してポリガンマ関数と呼びます．

　引数が整数のときのジガンマ・トリガンマ関数

φ（ｎ）＝Σ1/kｰγ，φ’（ｎ）＝π^2/6-Σ1/k^2

また，

φ^(k)（x）=(-1)^(k+1)k!Σ1/(n+x)^(k+1)より，ｘ＝１におけるポリガンマ関数値は，

φ(1)=-γ，φ'(1)=π2/6=ζ(2)，φ^(k)（1）=(-1)^(k+1)k!ζ(k+1)

になります．

（ガンマ関数・ポリガンマ関数のサブルーチン）

　Γ（ｎ＋１）＝ｎ！より，ガンマ関数は階乗の補間関数であり，初等的でない関数の中で最も簡単かつ重要な数学的関数といえます．

　また，引数の値が半整数のときには，Γ（ｎ＋１／２）＝√π・（２ｎ）！／２^2nｎ！です．なお，ガンマ関数Γ（ｘ）はｘ＞０について微分可能で，ｘ＝１．４６１６３２１４４９・・・で最小となります．

　　Γ(1)＝１，Γ(1/2)＝√π

であることを知っていればたいてい間に合いますが，ガンマ関数は結構いろいろなところで出てきますから，ぜひ計算ルーチンを用意しておきたいものです．引数の値が整数（ｎ＋１）でも半整数（ｎ＋１／２）でもないガンマ関数に対する正確な計算は複雑ですが，物理・化学・工学などの諸分野の問題の解を求めるときには，どうしてもそれらの解を数値を用いて表す必要がでてきます．現在では膨大な数表を用いることは稀であって，コンピュータが広く使用されていますから，必要に応じて計算する方法がとられています．

　ガンマ関数Γ（ｘ）の値を求めるには，いろいろな近似式があり，例えば，ベルヌーイ数Ｂnを含む漸近展開

　　ｌｏｇΓ（ｘ）＝ｘｌｏｇｘ−ｘ＋１／２ｌｏｇ２π／ｘ＋Σ（−１）n-1Ｂn/(2n)(2n-1)x2n-1

は物理学（量子論）との関わりで，しばしば用いられています．

　このほかにも，ヘイスティングス，コリンジ，ハートらによる多項式展開などがありますが，ポリガンマ関数とのつながりを考えると，ヘイスティングスの多項式展開が優れていると思われます．ここでは，ヘイスティングス（Hastings）の８次多項式近似を紹介することにします．

Γ（ｘ＋１）≒1-0.577191652x+0.988205891x^2-0.897056937x^3+0.918206857x^4-0.756704078x^5+0.482199394x^6-0.193527818x^7+0.035868343x^8 （岩波全書「数学公式�V」森口，宇田川，一松著）．

　この近似式は有効範囲が０≦ｘ≦１に限られていますが，ガンマ関数にはΓ（ｘ＋１）＝ｘΓ（ｘ）の関係があり，この漸化式を繰り返し適用して，０≦ｘ≦１の範囲になるように再帰的な方法を用いればｘ≧１にも拡張することができます．また，ガンマ関数Γ（ｘ）はｘが大きくなるとすぐにオーバーフローエラーを起こしてしまいますので，プログラムではｌｎ（Γ（ｘ））＝ｌｎ（Γ（ｘ＋１））−ｌｎ（ｘ）によってその対数ｌｎ（Γ（ｘ））を求めて飽和現象を回避しています．

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

(4)ベータ関数

　ベータ関数（オイラーの第１種積分）は，

　　B(a,b)=∫(0,1)t^(a-1)(1-t)^(b-1)dt t0-1

によって定義されます．

　ここで，積分変数をｔからu=(1-t)/tによってuに変えると，

　　B(a,b)=∫(0,∞)u^(a-1)/(1+u)^(a+b)du u 0−∞

が得られます．

　ベータ関数とガンマ関数との間には

　　B(a,b)=Γ(a)Γ(b)/Γ(a+b)

の関係がありますから，ベータ関数はガンマ関数の兄弟分にあたります．

　　Γ(1/2)＝√π

を得るにはベータ関数が用いられます．この関数において，t=sin^2θとおくと，dt=2sinθcosθdθですから

　　B(a,b)=∫(0,1)t^(a-1)(1-t)^(b-1)dt=2∫(0,π/2)sin^(2a-1)θcos^(2b-1)θdθ

　ここで，a=1/2,b=1/2とすると

　　B(1/2,1/2)=2∫(0,π/2)dθ=π

　　Γ^2(1/2)/Γ(1)=π

Γ(1)=1であるからΓ(1/2)＝√πとなります．

　上式を一般化すると，

　　∫(a-b)(x-a)^m(b-x)^ndx=m!n!/(m+n+1)!(b-a)^(m+n+1)

が得られます．これらは受験参考書に必ず書いてある

　　∫(a-b)(x-a)(x-b)dx=-1/6(b-a)^3

　　∫(a-b)(x-a)(x-b)^2　dx=1/12(b-a)^4

という公式の一般化になっています．

　また，ジガンマ関数・トリガンマ関数との関係で，公式集にも収録されていないものの，応用上とくに重要な積分公式を次に示しておきます．

1/Γ(a)*∫(0,∞)logt*t^(a-1)exp(-t)dt=φ(a)

1/Γ(a)*∫(0,∞)(logt)^2*t^(a-1)exp(-t)dt=φ(a)^2+φ'(a)

1/B(a,b)*∫(0,1)logt*t^(a-1)(1-t)^(b-1)dt=-φ(a)+φ(a+b)

1/B(a,b)*∫(0,1)(logt)^2*t^(a-1)(1-t)^(b-1)dt={-φ(a)+φ(a+b)}^2+φ'(a)-φ'(a+b)

1/B(a,b)*∫(0,∞)logt*t^(a-1)/(1+t)^(a+b)dt=φ(a)-φ(b)

1/B(a,b)*∫(0,∞)(logt)^2*t^(a-1)/(1+t)^(a+b)dt={φ(a)-φ(b)}^2+φ'(a)+φ'(b)

∫(0,∞)logt*t^(a-1)exp(-t)dt=Γ'(a)

∫(0,∞)(logt)^n*t^(a-1)exp(-t)dt=Γ(n)(a)

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝