確率論では，h(x)=exp(itx)とおいた場合の期待値をｔの関数φ(t)とみて，特性関数といいます．

　　φ(t)=E[exp(itx)]=∫(-∞,∞)f(t)exp(itx)dt

フーリエ変換のカーネルはh(x)=exp(-itx)ですから，特性関数は１種のフーリエ変換（＋ｉフーリエ変換）と考えることができます．→【補】

　

　特性関数はすべての確率分布に対して存在し，正規分布Ｎ（μ，σ^2）の特性関数は

　　exp(iμt-σ^2t^2/2)

コーシー分布の特性関数は

　　exp(iαt-|t|β)

と表されます．また，畳み込みのフーリエ変換はフーリエ変換の単なる積になりますから，畳み込みの特性関数はそれぞれの分布の特性関数の積

　　φx+y(t)=φx(t)*φy(t)

で表されます．→【補】

　

　ガウス関数（Gaussian）とローレンツ関数（Lorentzian）の畳み込みは，これを逆変換

　　1/2π∫(-∞,∞)φx+y(t)exp(-itx)dt

することによって，積分式

　　h(x)=1/π∫(0,∞)exp(-σ^2t^2/2-βt)cos{(γ-x)t}dt

で表すことができます．この関数はVoigt関数と呼ばれ，分光学の分野では非常に重要な関数となっています．

　

　しかし，フォークト関数の密度関数は積分関数を含んでおり，このタイプは，通常，解析的に積分の結果を得ることは困難です．かえって，積分表示のままで要求される精度の演算を行い，うまくいく場合もありますから，数値積分的に求めることでよい場合は，ガウス・エルミート法が適していると思われます．

　

　ガウス・エルミート則については，ご存じかとも思いますが，エルミート関数の零点を分点とし，分点におけるf(t)の値に重みを掛け，補間する公式による数値積分計算で，台形則やシンプソン則を拡張したものです．フォークト関数そのものの形は素直な関数ですから期待がもてそうですし，実際，ガウス・エルミート法によるVoigt関数の計算ソフトもあり，収束良く求まりますから，やはり，ガウス・エルミート法による数値積分がベストの選択と思われました．

　

　ところが，東北大学・薬学部の外山先生が所望しているのは，Voigt関数を数値的にではなくて，解析的に表現するのが目的であり，その希望に応えるべく，超幾何関数，放物柱関数，あるいは，複素誤差関数を用いた計算法を考えることにしました．

　

　そのすべてを一度にここに記すには無理があり，今回のコラムでは，超幾何関数による計算プログラムを紹介したいと思います．超幾何関数はその微分も積分も超幾何関数で表せるので非常に便利ですし，解析的な解を模索したいのであれば，超幾何関数法は有用と考えられました．

　

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

　

【１】超幾何関数を用いたVoigt関数の計算

　

　実誤差関数は，超幾何関数を用いた方法で級数展開することによって，精度よく計算できます．したがって，Voigt関数が複素誤差関数で表されると知ったとき，真っ先に浮かんだのが超幾何関数を用いる方法です．

　

　そこで，

　　f(t)=exp(-at^2-bt)

　　g(t,z)=cos{(c-z)t}

とおくとき,

　　h(z)=∫(0,∞)f(t)g(t,z)dt　　 （a,b,cは正の定数）

なる関数h(z)が，もし，合流型あるいはガウス形超幾何関数の形に表せたなら，超幾何関数には微分・積分してもふたたび超幾何関数になるという特性

　　d/dx1F1(a,c,x)=a/c1F1(a+1,c+1,x)

　　d/dx2F1(a,b,c,x)=ab/c2F1(a+1,b+1,c+1,x)

があり，微分・積分計算が非常に楽になりますから大変都合がよい．ところが，f(t),g(t)いずれも簡単な超幾何関数になるのものの，それらの積の広義積分については，公式集を探ってみても自力ではうまく超幾何関数に変形できませんでした．

　

　h(z)は超幾何級数展開できないのではと途方にくれていたところ，鈴鹿高専の奥井重彦先生より，新数学公式集（丸善）のp451,(2.5.36.2)式：

　　∫(0,∞)x^(α-1)exp(-ax^2)cos(cx)dx

　　=1/2a^(α/2)Γ(α/2)1F1(α/2,1/2;-c^2/4a)

が利用できるという貴重な情報が寄せられました．

　

　　exp(-bx)=Σ(ｰbx)^k/k!

と展開すれば，

　　∫(0,∞)exp(-ax^2-bx)cos(cx)dx

　　∫(0,∞)exp(-ax^2)cos(cx)Σ(ｰbx)^k/k!dx

　　=Σ(-b)^k/k!∫(0,∞)x^kexp(-ax^2)cos(cx)dx

　　=Σ(ｰb)^k/k!/2a^((k+1)/2)Γ((k+1)/2)1F1((k+1)/2,1/2;-c^2/4a)

となります．ここでは項別積分や順序交換に関するルベーグ積分の重要な定理が使われていることに注意して下さい．この超幾何関数を用いる式を利用すれば，うまく計算できそうです．

　

1000 '

1010 ' *** Voigt 関数の計算（超幾何関数） ***

1020 '

1030 *VOIGT.INIT:

1040 '[in]

1050 'x : [argument]

1060 'p1: [weight]

1070 'p2: [location]

1080 'p3: [scale of Gaussian]

1090 'p4: [scale of Lorentzian]

1100 '

1110 '[out]

1120 'voigt#: f(x,p2,p3,p4)

1130 'v1#,v2#,v3#,v4#: Jacobian of p1*f(x,p2,p3,p4)

1140 'vd#: df/dx

1150 'vi#: integral(-∞,x)f(t)dt

1160 '

1170 PI#=3.141592653589732#

1180 QI#=SQR(PI#)

1190 '

1200 Q1=P3*P3/2

1210 Q2=P2-X

1220 Q3=P4

1230 UUX=-Q2*Q2/Q1/4

1240 UUQ=ABS(Q2)

1250 RETURN

1260 '

1270 ' *** Voigt 関数 ***

1280 '

1290 *VOIGT.HGF:

1300 GOSUB *VOIGT.INIT

1310 Z=UUX

1320 J=-2:GOSUB *VOIGT.HGF.SUB:G1#=G#

1330 J=-1:　　　　　　　　　　:G2#=1#

1340 J=0 :GOSUB *VOIGT.HGF.SUB:G3#=G#

1350 J=1 :GOSUB *VOIGT.HGF.SUB:G4#=G#

1360 '

1370 EPSL=1E-09

1380 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

1390 W1#= QI#/SQR(Q1)/2

1400 W2#=-1/Q1/2*Q3

1410 SS#=0:S0#=W1#*H3#+W2#*H4#

1420 J=0

1430 WHILE ABS((SS#-S0#)/S0#)>EPSL

1440　 A=(1+J)/2:C=1/2

1450　 H5#=(C-A)/A*H1#+(A*2-C+Z)/A*H3#

1460　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*A

1470　 '

1480　 J=J+1

1490　 A=(1+J)/2:C=1/2

1500　 H6#=(C-A)/A*H2#+(A*2-C+Z)/A*H4#

1510　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*A

1520　 '

1530　 H1#=H3#:H3#=H5#

1540　 H2#=H4#:H4#=H6#

1550　 SS#=S0#

1560　 S0#=S0#+W1#*H5#+W2#*H6#

1570　 J=J+1

1580 WEND

1590 VOIGT#=S0#/PI#

1600 RETURN

1610 '

1620 *VOIGT.HGF.SUB:

1630 A=(1+J)/2:C=1/2

1640 GOSUB *CONFLUENT.HGF

1650 RETURN

　

3450 '

3460 ' *** 合流型超幾何関数 ***

3470 '

3480 *CONFLUENT.HGF:' [Kummer]

3490 EPSL=1E-09

3500 G0#=0:G#=1

3510 T#=1:K=1

3520 WHILE ABS(G0#-G#)>EPSL

3530　 T#=T#*A/C*Z/K

3540　 G0#=G#

3550　 G#=G#+T#

3560　 K=K+1

3570　 A=A+1

3580　 C=C+1

3590 WEND

3600 G=G#

3610 RETURN

　

　3450行〜3610行が超幾何関数の級数和を求めるサブルーチンです．PC-9801の場合，単精度計算は１０進法で７桁の精度で格納され，６桁までの桁数で印刷されます（倍精度計算は格納・表示とも１６桁）．そこでプログラムでは，要求する精度が1E-9以下になったら反復を終了させるようにしています（3490行）．

　

　計算全体の精度を決めるのが３４９０行のεです．超幾何関数の計算精度は無限級数を有限回で打ち切ることになりますし，打ち切り誤差の他に，丸め誤差も入ってきますから，倍精度で計算しても，単精度と倍精度の中間くらいになるはずです．そのことを考慮に入れて，必要な精度に定めて下さい．

　

　このように，超幾何関数を用いれば確率計算を任意の精度で行うことができるという精度管理上のメリットが得られます．また，必要な精度があらかじめ定まっている場合には，その精度に見合った値を選んで使うことになります．

　

　また，フォークト関数の計算は，

　　∫(0,∞)exp(-ax^2-bx)cos(cx)dx

　　=Σ(ｰb)^k/k!/2a^((k+1)/2)Γ((k+1)/2)1F1((k+1)/2,1/2;-c^2/4a)

　　=w1f1+w2f2+w3f3+･･･　　 （fiは超幾何関数,wiはその重み）

のように無限級数の形になるわけですが，超幾何関数が無限級数であるうえに，それらの無限加重和ですから，フォークト関数の分布関数は二重の無限級数の形になっています．∞^2回の加重が必要ですから，大層時間のかかるうえに，計算誤差が大きいことが予想され，大変始末が悪いと考えられます．

　

　しかし，この不安は払拭されます．1260行〜1650行がフォークト関数の計算ルーチンですが，ここではガンマ関数の漸化式と超幾何関数のパラメータに関する３項漸化式を用いることにします．

　

　ガンマ関数Γ（ｘ）には，

　　Γ（ｘ＋１）＝ｘΓ（ｘ）

の関係があり，次のような漸化式が成り立ちます．

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

　

　したがって，ｘが正の整数ｎのときには，Γ（ｎ＋１）＝ｎ！が成り立ち，ガンマ関数は階乗の一般形となっていることがわかります．階乗の解析的補間をしている関数がガンマ関数なのです．

　

　このことから，w(i-2)がわかれば，wiは任意のｉに対して

　　wi/w(i-2)=Γ((i+1)/2)/Γ((i-1)/2)=(i-1)/2

と漸化式の形で書けます（1460,1510行）．wi/w(i-1)の計算は面倒なので，奇数番目と偶数番目を分けて考えた方がよいでしょう．なお，その際，

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

であることを知っておく必要があります（1390,1400行）．

　

　また，1450,1500行が，超幾何関数のパラメータに関する３項漸化式

　　1F1(a+1,c,x)=(c-a)/a1F1(a-1,c,x)+(2a-c+x)/a1F1(a,c,x)

　→奥井重彦「特殊関数とその応用」（森北出版），p104,(25)式

です．この漸化式を駆使することによって，超幾何関数値計算ルーチンへのアクセスを最大でも４回までに抑える工夫をしています（1320〜1350行）．

　

　すなわち，

　　wifi/w(i-2)f(i-2)

の比は簡単な形になりますから，w1f1→w3f3とw2f2→w4f4の値の計算から始めて，あとは漸化式を利用して，

　　w1f1→w3f3→w5f5→・・・

　　w2f2→w4f4→w6f6→・・・

と芋づる式に計算を進めていけばよいことになります．このような比をとるアイディアのおかげで計算所要時間をかなり短縮することができました．

　

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

　

【２】ヤコビアンの計算

　

　フォークト関数を

　　f(x,a,b,c,d)=∫(0,∞)exp(-at^2-bt)cos((c-x)t)dt

　　=Σ(ｰb)^k/k!/2a^((k+1)/2)Γ((k+1)/2)1F1((k+1)/2,1/2;-c^2/4a)

とおいて，そのヤコビアン：

　　df/da,df/db,df/dc

を求めてみましょう．

　

　その際，通常行われるのは右辺：

　　Σ(ｰb)^k/k!/2a^((k+1)/2)Γ((k+1)/2)1F1((k+1)/2,1/2;-c^2/4a)

をパラメータで微分する方法ですが，それでは計算が過度に複雑になってしまいますし，発展性もありません．そこで，再びルベーグ積分を利用することにします．

　

　ルベーグ積分の定理でよく使われるものを，成書から拾い上げてみると，

　　(1)項別積分に関するルベーグの優収束定理

　　(2)積分の順序交換に関するフビニの定理

　　(3)微積分の基本定理の一般化であるルベーグの微分定理

　　(4)ラドン・ニコディムの定理

などが列挙されますが，たとえば，

　　h(z)=∫(0,∞)exp(-at^2-bt)cos{(c-z)t}dt

なる関数において，df/dzを求めるという問題を考えてみましょう．

　

　このとき，微分と積分の順序を入れ替えて，

　　f'(z)=d{∫(0,∞)exp(-at^2-bt)cos{(c-z)t}dt}/dz

　　=∫(0,∞)f(t){dcos{(c-z)t}/dz}dt

　　=∫(0,∞)f(t)sin{(c-z)t}tdt

が，ルベーグ積分の定理で保証されています．

　

　一般に，f(x,a)がxに関して可積分，適当な積分可能な関数ψ(x)>=0があって，

　　|df(x,a)/da| <= ψ(x)

ならば，df(x,a)/daは可積分であり，積分と微分の順序を交換してよいのですが，上の例の場合，

　　|f(t){dg(t,z)/dz|<=Const f(t)

なので，OKということになります．ルベーグ積分の専門書の証明は難解ですが，得られた定理は実用的であり，応用数学で無意識のうちに使われています．

　

　これにより，ヤコビアンの計算はフォークト関数の計算のパラメータを変えるだけになり，1260行〜1650行の数カ所を手直しするだけで簡単に作成することができます．

　

1660 '

1670 ' *** ヤコビアン ***

1680 '

1690 *VOIGT.JACOBI1:

1700 GOSUB *VOIGT.HGF

1710 J1#=VOIGT#

1720 RETURN

1730 '

1740 ' *** ヤコビアン ***

1750 '

1760 *VOIGT.JACOBI2:

1770 GOSUB *VOIGT.INIT

1780 IF Q2=0 THEN J2#=0#:RETURN

1790 Z=UUX

1800 J=-2:GOSUB *VOIGT.JACOBI2.SUB:G1#=G#

1810 J=-1:GOSUB *VOIGT.JACOBI2.SUB:G2#=G#

1820 J=0 :GOSUB *VOIGT.JACOBI2.SUB:G3#=G#

1830 J=1 :GOSUB *VOIGT.JACOBI2.SUB:G4#=G#

1840 '

1850 EPSL=1E-09

1860 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

1870 W1#=-QI#/2/Q1^(3/2)*Q2/2

1880 W2#= 1/Q1^2*Q2/2*Q3

1890 SS#=0:S0#=W1#*H3#+W2#*H4#

1900 J=0

1910 WHILE ABS((SS#-S0#)/S0#)>EPSL

1920　 A=(3+J)/2:C=3/2

1930　 H5#=(C-A)/A*H1#+(A*2-C+Z)/A*H3#

1940　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*A

1950　 '

1960　 J=J+1

1970　 A=(3+J)/2:C=3/2

1980　 H6#=(C-A)/A*H2#+(A*2-C+Z)/A*H4#

1990　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*A

2000　 '

2010　 H1#=H3#:H3#=H5#

2020　 H2#=H4#:H4#=H6#

2030　 SS#=S0#

2040　 S0#=S0#+W1#*H5#+W2#*H6#

2050　 J=J+1

2060 WEND

2070 J2#=S0#/PI#*P1

2080 RETURN

2090 '

2100 *VOIGT.JACOBI2.SUB:

2110 A=(3+J)/2:C=3/2

2120 GOSUB *CONFLUENT.HGF

2130 RETURN

2140 '

2150 ' *** ヤコビアン ***

2160 '

2170 *VOIGT.JACOBI3:

2180 GOSUB *VOIGT.INIT

2190 Z=UUX

2200 J=-2:GOSUB *VOIGT.JACOBI3.SUB:G1#=G#

2210 J=-1:GOSUB *VOIGT.JACOBI3.SUB:G2#=G#

2220 J=0 :GOSUB *VOIGT.JACOBI3.SUB:G3#=G#

2230 J=1 :GOSUB *VOIGT.JACOBI3.SUB:G4#=G#

2240 '

2250 EPSL=1E-09

2260 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

2270 W1#=-QI#/2/Q1^(3/2)/2*P3

2280 W2#= 1/Q1^2/2*Q3*P3

2290 SS#=0:S0#=W1#*H3#+W2#*H4#

2300 J=0

2310 WHILE ABS((SS#-S0#)/S0#)>EPSL

2320　 A=(3+J)/2:C=1/2

2330　 H5#=(C-A)/A*H1#+(A*2-C+Z)/A*H3#

2340　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*A

2350　 '

2360　 J=J+1

2370　 A=(3+J)/2:C=1/2

2380　 H6#=(C-A)/A*H2#+(A*2-C+Z)/A*H4#

2390　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*A

2400　 '

2410　 H1#=H3#:H3#=H5#

2420　 H2#=H4#:H4#=H6#

2430　 SS#=S0#

2440　 S0#=S0#+W1#*H5#+W2#*H6#

2450　 J=J+1

2460 WEND

2470 J3#=S0#/PI#*P1

2480 RETURN

2490 '

2500 *VOIGT.JACOBI3.SUB:

2510 A=(3+J)/2:C=1/2

2520 GOSUB *CONFLUENT.HGF

2530 RETURN

2540 '

2550 ' *** ヤコビアン ***

2560 '

2570 *VOIGT.JACOBI4:

2580 GOSUB *VOIGT.INIT

2590 Z=UUX

2600 J=-2:　　　　　　　　　　　　:G1#=1#

2610 J=-1:GOSUB *VOIGT.JACOBI4.SUB:G2#=G#

2620 J=0 :GOSUB *VOIGT.JACOBI4.SUB:G3#=G#

2630 J=1 :GOSUB *VOIGT.JACOBI4.SUB:G4#=G#

2640 '

2650 EPSL=1E-09

2660 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

2670 W1#=-1/Q1/2

2680 W2#= QI#/2/Q1^(3/2)/2*Q3

2690 SS#=0:S0#=W1#*H3#+W2#*H4#

2700 J=0

2710 WHILE ABS((SS#-S0#)/S0#)>EPSL

2720　 A=(2+J)/2:C=1/2

2730　 H5#=(C-A)/A*H1#+(A*2-C+Z)/A*H3#

2740　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*A

2750　 '

2760　 J=J+1

2770　 A=(2+J)/2:C=1/2

2780　 H6#=(C-A)/A*H2#+(A*2-C+Z)/A*H4#

2790　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*A

2800　 '

2810　 H1#=H3#:H3#=H5#

2820　 H2#=H4#:H4#=H6#

2830　 SS#=S0#

2840　 S0#=S0#+W1#*H5#+W2#*H6#

2850　 J=J+1

2860 WEND

2870 J4#=S0#/PI#*P1

2880 RETURN

2890 '

2900 *VOIGT.JACOBI4.SUB:

2910 A=(2+J)/2:C=1/2

2920 GOSUB *CONFLUENT.HGF

2930 RETURN

2940 '

2950 ' *** Voigt 関数の微分 ***

2960 '

2970 *VOIGT.DIFFER:

2980 GOSUB *VOIGT.INIT

2990 IF Q2=0 THEN VD#=0#:RETURN

3000 GOSUB *VOIGT.JACOBI2

3010 VD#=-S0#/PI#

3020 RETURN

3030 '

3040 ' *** Voigt 関数の積分 ***

3050 '

3060 *VOIGT.INTEGR:

3070 GOSUB *VOIGT.INIT

3080 IF Q2=0 THEN VI#=.5#:RETURN

3090 Z=UUX

3100 J=-2:GOSUB *VOIGT.INTEGR.SUB:G1#=G#

3110 J=-1:　　　　　　　　　　　 :G2#=1#

3120 J=0 :GOSUB *VOIGT.INTEGR.SUB:G3#=G#

3130 J=1 :GOSUB *VOIGT.INTEGR.SUB:G4#=G#

3140 '

3150 EPSL=1E-09

3160 H1#=G1#:H2#=G2#:H3#=G3#:H4#=G4#

3170 W1#= QI#/SQR(Q1)/2*UUQ

3180 W2#=-1/Q1/2*Q3*UUQ

3190 SS#=0:S0#=W1#*H3#+W2#*H4#

3200 J=0

3210 WHILE ABS((SS#-S0#)/S0#)>EPSL

3220　 A=(1+J)/2:C=3/2

3230　 H5#=(C-A)/A*H1#+(A*2-C+Z)/A*H3#

3240　 W1#=W1#*Q3*Q3/(J+1)/(J+2)/Q1*A

3250　 '

3260　 J=J+1

3270　 A=(1+J)/2:C=3/2

3280　 H6#=(C-A)/A*H2#+(A*2-C+Z)/A*H4#

3290　 W2#=W2#*Q3*Q3/(J+1)/(J+2)/Q1*A

3300　 '

3310　 H1#=H3#:H3#=H5#

3320　 H2#=H4#:H4#=H6#

3330　 SS#=S0#

3340　 S0#=S0#+W1#*H5#+W2#*H6#

3350　 J=J+1

3360 WEND

3370 V#=S0#/PI#

3380 VI#=.5#-V#*SGN(Q2)

3390 RETURN

3400 '

3410 *VOIGT.INTEGR.SUB:

3420 A=(1+J)/2:C=3/2

3430 GOSUB *CONFLUENT.HGF

3440 RETURN

　

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

　

【３】まとめ

　

　超幾何関数を用いたVoigt関数の計算プログラムを作成しました．この方法はシンプソン則やガウス・エルミート則などの数値積分法とは違って，（打ち切り誤差と丸め誤差を除いて）解析的な解が求まります．

　

　超幾何関数が無限級数である上に，超幾何関数の無限関数列になっているのですが，超幾何関数に最大でも４回アクセスするだけで済みますから，速度的にも満足できるものでした．あとは要求された精度を満たすようにεの値を決めてください．

　

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

　

【補】積率母関数と特性関数

　

　関数f(x)としてsinxやxの無限積分を考えると前者は不定，後者は発散してしまいます．一方，これらの関数にexp(-x)を掛け合わせたf(x)exp(-x)を無限積分すると収束します．x>0のとき，ガンマ関数

　　Γ(s)=∫(0-∞）x^(s-1)exp(-x)dx

の存在が知られているルーツにもこのような理由があるからです．

　

　すなわち，exp(-x)は無限積分において不定や発散する関数を収束させる働きをもっていることが理解されます．このことより，exp(-x)の代わりにもうひとつの変数ｓを含んだexp(-sx)を考え，

　　F(s)=∫(-∞,∞）f(x)exp(-sx)dx

とおくと，無限積分の後，ｘの関数はｓの関数に変換されます．この操作をラプラス変換と呼びます．ラプラス変換において変数ｓは複素変数であり，フーリエ変換はラプラス変換におけるパラメータｓの実部が０である場合に相当します．

　

［１］積率母関数

　

　h(x)=exp(tx)とおいた場合の期待値をｔの関数M(t)とみて，積率母関数といいます．

　　M(t)=E[exp(tx)]=∫(-∞,∞)f(t)exp(tx)dt

ラプラス変換のカーネルはh(x)=exp(-tx)ですから，積率母関数は１種のラプラス変換（＋ラプラス変換）と考えることができます．

　

　また，

　　M(t)=E[etx]=E[1+tx+(tx)^2/2!+(tx)^3/3+･･･]

　　=1+tE[x]+t^2/2!E[x^2]+t^3/3!E[x^3]+･･･

より，

　　M(t)=Σμ'k(t)^k/k!　　　k=0,1,2,･･･

を得ます．すなわち，積率母関数は原点まわりの積率の指数型母関数になっていることが理解されます．積率母関数を微分してｔ＝０とおくことによって原点まわりの積率μ'kが順次求まります．

　　d(k)M/dt(k)|t=0=μ'k=E[x^k]

　

（例題）正規分布の積率母関数は，M(t)=exp(μt+σ2t2/2)より，μ1，μ2を求めよ．

　

　M'(t)=(μ+σ2t)exp(μt+σ2t2/2)より，E[x]=M'(0)=μ

　M"(t)=(σ2+(μ+σ2t)^2)exp(μt+σ2t2/2)より，E[x^2]=M"(0)=σ2+μ2

を得ることができます．したがって，μ1=μ,μ2=σ2

　

　積率母関数には，「和の分布の積率母関数は積率母関数の積で表される．」という重要な性質があります．すなわち，x1,x2,...,xnが独立で，それぞれの積率母関数をMx1(t),Mx2(t),･･･,Mxn(t)とするとy=x1+x2+･･･+xnの積率母関数My(t)は

　　My(t)=ΠMxi(t)

で表されるというものです．とくに，x1,x2,･･･,xnの積率母関数が同じを積率母関数Mx(t)をもつとき，My(t)=[Mx(t)]^nとなります．

　

　正規分布の和の分布について考えてみましょう．ｘがN(μx,σx2)に，ＹがN(μy,σy2)にしたがい，両者が独立であれば，ｘ＋ｙの積率母関数は

　　Mx+y(t)=Mx(t)*My(t)=exp(μxt+σx2t2/2)exp(μyt+σy2t2/2)

　　=exp((μx+μy)t+(σx2+σy2)t2/2)

　

　これはＮ（μx+μy，σx2+σy2）の積率母関数にほかなりません．したがって，正規分布の和の分布はまた正規分布となります．これを正規分布の再生性といいます．なお，ポアソン分布や負の２項分布，コーシー分布やガンマ分布も再生性を有しています．

　

　一方，差の分布の積率母関数は，

　　Mx-y(t)=Mx(t)*My(-t)

で表されます．例題と同様に，正規分布の差の分布は

　　Mx-y(t)=Mx(t)*My(-t)=exp(μxt+σx2t2/2)exp(-μyt+σy2t2/2)

　　=exp((μx-μy)t+(σx2+σy2)t2/2)

すなわち，Ｎ（μx-μy，σx2+σy2）の正規分布になることを示すことができます．ところが，ポアソン分布の差の分布はポアソン分布にはならず，ベッセル関数を用いて表される形になります．

　

［２］特性関数

　

　積率母関数の欠点は，積率をもたないコーシー分布やブラウンノイズ関数などに対しては積率母関数が定義されないということです．そこで，複素関数を導入して

　　h(x)=exp(itx)

としたものが特性関数です．

　　φ(t)=E[exp(itx)]=∫(-∞,∞)f(t)exp(itx)dt

　

　フーリエ変換のカーネルはh(x)=exp(-itx)ですから，特性関数は1種のフーリエ変換（＋ｉフーリエ変換）と考えることができます．また，上式において，exp(itx)にオイラーの公式

　　exp(itx)=cos(tx)+isin(tx)

を適用すると，特性関数は次のように表現できます．

　　φ(t)=∫(-∞,∞)f(t)cos(tx)dt+i∫(-∞,∞)f(t)sin(tx)dt

　

　正規分布Ｎ（μ，σ2）の積率母関数は，M(t)=exp(μt+σ2t2/2)ですから，その特性関数はexp(iμt-σ2t2/2)となります．また，特性関数はすべての確率分布に対して存在し，コーシー分布の特性関数はexp(iμt-|t|σと表されます．

　

　また，畳み込みのフーリエ変換はフーリエ変換の単なる積になりますから，畳込みの特性関数はそれぞれの分布の特性関数の積

　　φx+y(t)=φx(t)*φy(t)

で表されます．また，差の分布の特性関数は，

　　φx-y(t)=φx(t)*φy(-t)

で表されます．

　

　積率に関しても，積率母関数と同様なことは特性関数でもいえて，

　　φ(t)=Σμ'r(it)^r/r!

　　μ'r=(-i)^rd(r)φ/dt(r)|t=0

が示されます．

　

（例題）特性関数も積率母関数同様に和の分布を求めるときなど利用されていますが，ここでは，区間(0,1)の一様分布の特性関数が

　　φ(t)=exp(it/2)sin(t/2)/(t/2)

となることを利用して，一様乱数ｒiをｎ個合計したものの分布が，ｎ→∞の極限で正規分布になることを示してみましょう．

　

　一様乱数をｎ個の合計のしたものの分布の特性関数は

　　[φ(t)]^n=exp(int/2){sin(t/2)/(t/2)}^n

　

一方，シンク関数

　　sinx/x=1-1/3!x2+5!x4-･･･（ベキ級数表示）

　　=Σ(-1)^mx^2m/(2m+1)!

の解が±π，±２π，±３π，±４π，・・・となることを利用して，無限積表示すると

　　sinx/x=(1-x2/π2)(1-x2/4π2)(1-x2/9π2)(1-x2/16π2)･･･

　　=Π(1-x2/k2π2)　　（無限積表示）

　

　ｋが大きいとき，

　　(1-x2/k2π2)〜exp(-x2/k2π2)

　　Π(1-x2/k2π2)〜exp(-x2/π2Σ1/k2)

ここで，Σ1/k2=π2/6=ζ(2)より，

　　Π(1-x2/k2π2)〜exp(-x2/6)

これより，{sin(t/2)/(t/2)}^n→exp(-nt2/24)ですから，

　　[φ(t)]^n→exp(int/2-nt2/24)

　

　正規分布Ｎ（μ，σ2）の特性関数はexp(iμt-σ2t2/2)ですから，この結果はｎ個の独立した一様乱数の和の分布は平均値ｎ／２，分散ｎ／１２の正規分布に近づくことを示していて，これは，一様分布の場合について中心極限定理を示したものとなっています．

　

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