今回のコラムでは，放物柱関数を用いたVoigt関数の計算法について紹介します．実は放物柱関数は超幾何関数に帰着するので，前回，前々回のプログラムと足したり引いたりする順番が異なるだけであって，本質的な違いはありません．そこで冒頭では，畳み込み（convolution）について，別の角度からながめてみたいと思います．

　

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

　

【１】分布の混合と畳み込み

　

　そのためには，まず，分布の混合（compoundまたはmixture of probability distribution）という用語について説明します．「分布の混合」は２つの意味で使われています．

　

　１つは，複数の独立な確率分布

　　f1(x),f2(x),･･･,fk(x)

がある比率

　　w1,w2,･･･,wk　　（Σwi=1，wi>0）

をもって混ぜ合わさっている場合であり，

　　f(x)=w1f1(x)+w2f2(x)+･･･+wkfk(x)=Σwifi(x)

が求める確率密度関数となります．これはいわば単純な加重和であって，混合される各分布は成分ということになります．また，合計１になる正の重み自体が１つの確率分布（g）に従うと考えることができます．例をあげると，非心χ2分布はχ2分布（f）とポアソン分布（g）の混合となっています．

　

　もう１つは，確率密度関数のパラメータが確率変数となっている確率分布関数を指す場合です．母数θを含む確率密度関数f(x)をパラメータまで含めてf(x,θ)と表します．f(x,θ)の母数θがまた確率密度関数g(θ)に従うとすると

　　h(x)=∫f(x,θ)g(θ)dθ

によって新しい確率密度関数h(x)が定義されます．また，混合分布の累積分布関数は

　　H(x)=∫h(x)dx=∫F(x,θ)g(θ)dθ=∫F(x,θ)dG(θ)

となります．

　

　これでθに関する重み付き平均を考えたことになりますが，位置母数モデルf(x,θ)=f(x-θ)のとき，hはfとgとの畳み込みそのものです．

　

　Ｘ線，γ線などの電磁波はそれぞれの線スペクトルに固有の幅と分布をもっていて，光の線スペクトルのようなコーシー分布（ｆ）：

　　f(x)=1/π･β/(β^2+(x-θ)^2)

を示すものを分光器で測定したとすると，分光器には固有の分解能があり，それは正規分布（ｇ）

　　g(θ)=1/√2πσexp(-θ^2/2σ^2)

で近似できることが多いわけですから，測定したスペクトルの分布ｈはコーシー分布と正規分布を合成したものになります．すなわち，フォークト関数（ｈ）はコーシー分布と正規分布の位置混合になっているというわけです．

　

　また，ある生体システム（ｇ）にあるパルス（ｆ）を入力した場合，出力としてｈ（ｔ）が観察された場合，たとえば，薬物経口投与後の血中濃度ｈ（ｔ）は，静注後の血中濃度ｆと吸収速度ｇの関数として表現することができます．このようにｆとｇを合成した関数は畳み込みと呼ばれ，一般に，分布のたたみ込みは

　　p*q(x)=∫(-∞,∞)p(x-t)q(t)dt

として定義されます．この式は，薬物体内動態評価のみならず，量子力学，数理生物学，赤外線分析，Ｘ線回折，クロマトや反応工学および多くの物理現象にわたって現れる普遍的関係式です．

　

　一方，尺度母数モデルf(x,θ)=f(x/θ)/θのとき，hはfとgの尺度混合であるといいます．g(θ)としてはガンマ分布（カイ２乗分布）が用いられることが多く，たとえば，ポアソン分布（f）をガンマ分布（g）で混合すれば，負の２項分布（h），指数分布を（f）をガンマ分布（g）で混合すれば，パレート分布（h）が得られます．ｔ分布も正規分布（f）のσ2がガンマ分布（g）にしたがうと仮定，すなわちガンマ分布で尺度混合した混合分布です．

　

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

　

【２】放物柱関数を用いたVoigt関数の計算

　

　超幾何関数を用いたVoigt関数の計算では，exp(-bt)を展開して，超幾何関数列としましたが，今回のコラムでは，cos{(c-x)t}を展開して,放物柱関数経由でVoigt関数の計算を行いました．

　

　用いた公式は，新数学公式集（丸善）のp343,(2.3.15.3)式：

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

　　=Γ(α/2)(2a)^(-α/2)exp(b^2/8a)D(-α,b/√2a)

に，ｃｏｓの無限ベキ級数表示

　　cosx=Σx^k/(2k+1)!

の展開式を組み合わせたものです．

　

　すなわち，

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

　　=Σ(c-x)^k/(2k+1)!∫(0,∞)t^kexp(-ax^2-bt)dt

　　=Σ(c-x)^k/(2k+1)!Γ((k+1)/2)(2a)^(-(k+1)/2)exp(b^2/8a)D(-k-1,b/√2a)

です．

　

　この式を利用すると，フォークト関数の積分値の計算が非常に簡単に行えるのですが，収束が若干不安定になりました．

　

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 Q4=EXP(-UUX)

1260 '

1270 UUX=P4*P4/P3/P3/2

1280 Q=Q2/P3

1290 QQ=Q*Q

1300 RETURN

1310 '

1320 ' *** Voigt 関数 ***

1330 '

1340 *VOIGT.HGF:

1350 GOSUB *VOIGT.INIT

1360 Z=UUX

1370 A=-1/2:C=1/2:GOSUB *CONFLUENT.HGF:G1#=G#

1380 A=0　 :C=3/2:　　　　　　　　　　:G2#=1#

1390 A=1/2 :C=1/2:GOSUB *CONFLUENT.HGF:G3#=G#

1400 A=1　 :C=3/2:GOSUB *CONFLUENT.HGF:G4#=G#

1410 '

1420 EPSL=1E-09

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

1440 W1#= QI#/SQR(2)/P3

1450 W2#=-P4/P3/P3

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

1470 J=0

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

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

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

1510　 W1#=-W1#*QQ/(J+1)/2

1520　 '

1530　 A=J+1　:C=3/2

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

1550　 W2#=-W2#*QQ/(J+1/2)/2

1560　 '

1570　 H1#=H3#:H3#=H5#

1580　 H2#=H4#:H4#=H6#

1590　 SS#=S0#

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

1610　 J=J+1

1620 WEND

1630 VOIGT#=S0#/PI#

1640 RETURN

1650 '

1660 ' *** ヤコビアン ***

1670 '

1680 *VOIGT.JACOBI1:

1690 GOSUB *VOIGT.HGF

1700 J1#=VOIGT#

1710 RETURN

1720 '

1730 ' *** ヤコビアン ***

1740 '

1750 *VOIGT.JACOBI2:

1760 GOSUB *VOIGT.INIT

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

1780 Z=UUX

1790 A=1/2 :C=1/2:GOSUB *CONFLUENT.HGF:G1#=G#

1800 A=1　 :C=3/2:GOSUB *CONFLUENT.HGF:G2#=G#

1810 A=3/2 :C=1/2:GOSUB *CONFLUENT.HGF:G3#=G#

1820 A=2　 :C=3/2:GOSUB *CONFLUENT.HGF:G4#=G#

1830 '

1840 EPSL=1E-09

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

1860 W1#=-QI#/SQR(2)/P3/P3*Q

1870 W2#= P4*2/P3/P3/P3*Q

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

1890 J=1

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

1910　 A=J+1/2:C=1/2

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

1930　 W1#=-W1#*QQ*(J+1)/J/(J+1)/2

1940　 '

1950　 A=J+1　:C=3/2

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

1970　 W2#=-W2#*QQ*(J+1)/J/(J+1/2)/2

1980　 '

1990　 H1#=H3#:H3#=H5#

2000　 H2#=H4#:H4#=H6#

2010　 SS#=S0#

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

2030　 J=J+1

2040 WEND

2050 J2#=S0#/PI#*P1

2060 RETURN

2070 '

2080 ' *** ヤコビアン ***

2090 '

2100 *VOIGT.JACOBI3:

2110 GOSUB *VOIGT.INIT

2120 Z=UUX

2130 A=1/2 :C=1/2:GOSUB *CONFLUENT.HGF:G1#=G#

2140 A=1　 :C=3/2:GOSUB *CONFLUENT.HGF:G2#=G#

2150 A=3/2 :C=1/2:GOSUB *CONFLUENT.HGF:G3#=G#

2160 A=2　 :C=3/2:GOSUB *CONFLUENT.HGF:G4#=G#

2170 '

2180 EPSL=1E-09

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

2200 W1#=-QI#/SQR(2)/P3/P3

2210 W2#= P4/P3/P3/P3*2

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

2230 J=0

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

2250　 A=J+3/2:C=1/2

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

2270　 W1#=-W1#*QQ*(J+2)/(J+1)*(J+3/2)/(J+1/2)/(J+2)/2

2280　 '

2290　 A=J+2　:C=3/2

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

2310　 W2#=-W2#*QQ*(J+2)/(J+1)*(J+3/2)/(J+1/2)/(J+3/2)/2

2320　 '

2330　 H1#=H3#:H3#=H5#

2340　 H2#=H4#:H4#=H6#

2350　 SS#=S0#

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

2370　 J=J+1

2380 WEND

2390 J3#=S0#/PI#*P1

2400 RETURN

2410 '

2420 ' *** ヤコビアン ***

2430 '

2440 *VOIGT.JACOBI4:

2450 GOSUB *VOIGT.INIT

2460 Z=UUX

2470 A=0　 :C=1/2:　　　　　　　　　　:G1#=1#

2480 A=1/2 :C=3/2:GOSUB *CONFLUENT.HGF:G2#=G#

2490 A=1　 :C=1/2:GOSUB *CONFLUENT.HGF:G3#=G#

2500 A=3/2 :C=3/2:GOSUB *CONFLUENT.HGF:G4#=G#

2510 '

2520 EPSL=1E-09

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

2540 W1#=-1/P3/P3

2550 W2#= QI#/SQR(2)*P4/P3/P3/P3

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

2570 J=0

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

2590　 A=J+1　:C=1/2

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

2610　 W1#=-W1#*QQ*(J+3/2)/(J+1/2)/(J+3/2)/2

2620　 '

2630　 A=J+3/2:C=3/2

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

2650　 W2#=-W2#*QQ*(J+3/2)/(J+1/2)/(J+1)/2

2660　 '

2670　 H1#=H3#:H3#=H5#

2680　 H2#=H4#:H4#=H6#

2690　 SS#=S0#

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

2710　 J=J+1

2720 WEND

2730 J4#=S0#/PI#*P1

2740 RETURN

2750 '

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

2770 '

2780 *VOIGT.DIFFER:

2790 GOSUB *VOIGT.INIT

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

2810 GOSUB *VOIGT.JACOBI2

2820 VD#=-S0#/PI#

2830 RETURN

2840 '

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

2860 '

2870 *VOIGT.INTEGR:

2880 GOSUB *VOIGT.INIT

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

2900 Z=UUX

2910 A=-1/2:C=1/2:GOSUB *CONFLUENT.HGF:G1#=G#

2920 A=0　 :C=3/2:　　　　　　　　　　:G2#=1#

2930 A=1/2 :C=1/2:GOSUB *CONFLUENT.HGF:G3#=G#

2940 A=1　 :C=3/2:GOSUB *CONFLUENT.HGF:G4#=G#

2950 '

2960 EPSL=1E-09

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

2980 W1#= QI#/SQR(2)/P3*UUQ

2990 W2#=-P4/P3/P3*UUQ

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

3010 J=0

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

3030　 A=J+1/2:C=1/2

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

3050　 W1#=-W1#*QQ*(J+1/2)/(J+3/2)/(J+1)/2

3060　 '

3070　 A=J+1　:C=3/2

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

3090　 W2#=-W2#*QQ*(J+1/2)/(J+3/2)/(J+1/2)/2

3100　 '

3110　 H1#=H3#:H3#=H5#

3120　 H2#=H4#:H4#=H6#

3130　 SS#=S0#

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

3150　 J=J+1

3160 WEND

3170 V#=S0#/PI#

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

3190 RETURN

3200 '

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

3220 '

3230 *CONFLUENT.HGF:' [Kummer]

3240 EPSL=1E-09

3250 G0#=0:G#=1

3260 T#=1:K=1

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

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

3290　 G0#=G#

3300　 G#=G#+T#

3310　 K=K+1

3320　 A=A+1

3330　 C=C+1

3340 WEND

3350 G=G#

3360 RETURN

　

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