（その４）の最後で，「非心χ2分布・非心ｔ分布・非心Ｆ分布は，直接的に超幾何関数で表すことはできないが，それぞれχ2分布・ｔ分布・Ｆ分布の加重和として間接的には表せるので，ここで述べた方法は非心χ2分布・ｔ分布・Ｆ分布に対する計算においても有用であるはずである．有用性が判明され次第，ＨＰにアップロードしてみたい．」と結びました．

　

　それを受けて，今回のコラムでは，非心χ2分布・非心ｔ分布・非心Ｆ分布の計算法とその計算プログラムを紹介します．一筋縄ではいかないことはハナから承知のうえだったのですが，思いの外，大変な作業となりました．しかし，その甲斐あって，やっとお披露目まで漕ぎ着けました．

　

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

　

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

　

　これ以降は前者の場合について考えてみましょう．

　

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

　

　非心χ2分布の分布関数を

　　Ｆ（ｘ，ν，λ）　　ν：自由度　　λ：非心パラメータ

とすると，

　　Ｆ（ｘ，ν，λ）＝Σexp(-λ/2)(λ/2)^j／Γ(j+1)Ｆ（ｘ，ν＋２ｊ，０）

　　　ｊ＝０〜∞

と表すことができます．

　

　Ｆ（ｘ，ν＋２ｊ，０）の部分は自由度ν＋２ｊのχ2分布の分布関数，exp(-λ/2)(λ/2)^j／Γ(j+1)はパラメータλ／２のポアソン確率を表しますから，Ｆ（ｘ，ν，λ）は，パラメータλ／２のポアソン確率をウェイトとするχ2分布の分布関数Ｆ（ｘ，ν＋２ｊ，０）の加重和に他なりません．これが，「非心χ2分布はχ2分布とポアソン分布の混合である」の意味です．

　

　同様に，非心Ｆ分布はＦ分布とポアソン分布の混合ですが，非心ｔ分布はｔ分布とポアソン分布の単純な混合にはなりません．ここでは詳しくは説明できませんが，ｔ分布自体が正規分布をガンマ分布で混合した混合正規分布であり，非心ｔ分布は，ｔ分布とポアソン分布，ｔ分布と離散正規分布（小生の造語）の複雑な混合になっているとだけ述べておきます．

　

　参考文献としては，応用統計の教科書として定評のある Johnson, Kotz らの「continuous univariate distributions-I,II」(John Wiley & Sons Inc.)をあげますが，それだけでは解らないでしょうから，以下にそのヒントを略記しました．なお，ｔ分布は左右対称でしたが，非心ｔ分布は対象形にはならないので，うまい手が使えないという厄介さもあります．

　

（ヒント）

　Σexp(-λ/2)(λ/2)^j／Γ(j+1)において，Σ(λ/2)^j／Γ(j+1)はexp(λ/2)のテイラー展開ですから，

　　Σexp(-λ/2)(λ/2)^j／Γ(j+1)=1

ここで，ｊにゲタをはかせて

　　Σexp(-λ/2)(λ/2)^(j+1/2)／Γ(j+3/2)

を考えてみます．直感的に

　　Σexp(-λ/2)(λ/2)^(j+1/2)／Γ(j+3/2)<1

と予想されますが，この式は解析的な形に表せるでしょうか？

　

　実は，上記の式は誤差関数の級数展開そのものであって，誤差関数の定義式を０の回りで級数展開し，項別積分することにより得られます．Mathematica,Mapleなどを用いれば自動的に解が得られますが，私のやり方はΣ(λ/2)^(j+1/2)／Γ(j+3/2)を超幾何関数の形で，

　　1Ｆ1(1,3/2,x)

に還元しておいてから，誤差関数Φ（ｘ）にたどり着くというものでした．この超幾何関数は，級数Σ(λ/2)^(j+1/2)／Γ(j+3/2)の項比が

　　ａj+1／ａj＝(j+1)/(j+3/2)･(λ/2)/(j+1)

であることから導出できます．

　

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

　

　前節で示したことより，早速，ひとつの問題が提起されます．χ2分布の分布関数は

　　Ｆ（ｘ）＝(x/2)^(ν/2)exp(-x/2)1Ｆ1(1,1+ν/2,x/2)

そのもの，一方，ｔ分布やＦ分布の分布関数はベータ分布と関係していて，

　　Ｆ（ｘ）＝１／ｐｘ^p（１−ｘ）^q2Ｆ1(p+q,1,p+1,x)／Β（ｐ，ｑ）

と表されます．

　

　いずれの分布関数も超幾何関数を含んでいますが，超幾何関数が無限級数

　　Ｆ（ｘ）＝Σａiｘ^i　　ｉ＝０〜∞

であるうえに，非心分布はそれらの分布の無限加重和

　　Ｇ（ｘ）＝ΣｗjＦj（ｘ）　　ｊ＝０〜∞

ですから，非心分布の分布関数は二重の無限級数の形

　　Ｇ（ｘ）＝ΣｗjΣａijｘ^iＦj（ｘ）　　ｉ＝０〜∞，ｊ＝０〜∞

になっています．

　

　前回掲げたプログラムでは，超幾何関数のパラメータに関する漸化式を駆使することによって，超幾何関数値計算ルーチンへのアクセスを１回に減らすことができることを紹介しましたが，χ2分布・ｔ分布・Ｆ分布においては，たとえ１回のアクセスであっても，非心分布になると無限回の加重が必要ですから大層時間のかかるうえに，計算誤差が大きいことが予想され，大変始末が悪いと考えられることです．

　

　そこで，今回掲げたプログラムでは，非心分布でも超幾何関数値計算ルーチンへのアクセスをたった１回に抑える工夫をしました．たとえば，ポアソン分布：ｐ（ｘ）＝exp(-λ)λ^x／ｘ！では，

　　ｐ（ｘ）＝exp(-λ)λ^(x-1)／（ｘ−１）！・λ／ｘ

　　　　　　＝ｐ（ｘ−１）・λ／ｘ

したがって，ｐ（ｘ−１）が計算済みであれば，この漸化式を利用してｐ（ｘ）が計算できるわけで，ｐ（０）＝exp(-λ)から始めて，ｐ（１），ｐ（２），・・・と次々に計算して求めていけばよいことになります．

　

　超幾何関数の計算についても，このような比をとる方法を用いて省エネ計算します．実際の計算では荷重となるポアソン分布の他に，超幾何関数の項比が加わり漸化式はかなり複雑となりますが，このアイディアのおかげで計算時間をかなり節約することができました．

　

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

　

【プログラム】

　

1000 '

1010 ' ** NORMAL DISTRIBUTION **

1020 '

1030 *NORMAL.PERCENT.HGF:

1040 PP=1-P0:GOSUB *NORMAL.PERCENT:UU=UUN

1050 PROB$="N":GOSUB *NEWTON.HGF

1060 UUN=X

1070 RETURN

1080 '

1090 ' ** CHI-SQUARE DISTRIBUTION **

1100 '

1110 *CHI2.PERCENT.HGF:

1120 PP=1-P0:GOSUB *CHI2.PERCENT:UU=UUX

1130 PROB$="X":GOSUB *NEWTON.HGF

1140 UUX=X

1150 RETURN

1160 '

1170 ' ** T-DISTRIBUTION **

1180 '

1190 *T.PERCENT.HGF:

1200 PP=1-P0:GOSUB *T.PERCENT:UU=UUT

1210 PROB$="T":GOSUB *NEWTON.HGF

1220 UUT=X

1230 RETURN

1240 '

1250 ' ** F-DISTRIBUTION **

1260 '

1270 *F.PERCENT.HGF:

1280 PP=1-P0:GOSUB *F.PERCENT:UU=UUF

1290 PROB$="F":GOSUB *NEWTON.HGF

1300 UUF=X

1310 RETURN

1320 '

1330 ' ** NON-CENTRAL CHI-SQUARE DISTRIBUTION **

1340 '

1350 *NCHI2.PERCENT.HGF:

1360 PP=1-P0:GOSUB *NCHI2.PERCENT1:UU=UUX

1370 PROB$="NX":GOSUB *NEWTON.HGF

1380 UUX=X

1390 RETURN

1400 '

1410 ' ** NON-CENTRAL T-DISTRIBUTION **

1420 '

1430 *NT.PERCENT.HGF:

1440 PP=1-P0:GOSUB *NT.PERCENT1:UU=UUT

1450 PROB$="NT":GOSUB *NEWTON.HGF

1460 UUT=X

1470 RETURN

1480 '

1490 ' ** NON-CENTRAL F-DISTRIBUTION **

1500 '

1510 *NF.PERCENT.HGF:

1520 PP=1-P0:GOSUB *NF.PERCENT1:UU=UUF

1530 PROB$="NF":GOSUB *NEWTON.HGF

1540 UUF=X

1550 RETURN

1560 '

1570 ' *** ニュートン法 ***

1580 '

1590 *NEWTON.HGF:

1600 EPSL=1E-09

1610 IMAX=50

1620 '

1630 Q0=P0

1640 D0=UU

1650 '

1660 SW=0

1670 ICNT=0

1680 WHILE SW=0

1690　 ICNT=ICNT+1:PRINT "*";

1700　 IF ICNT>IMAX THEN SW=1

1710　 DD=D0:GOSUB 1750:D0=DD

1720 WEND

1730 X=DD

1740 RETURN

1750 '

1760 ' ** 逐次計算 **

1770 '

1780 IF PROB$="N" THEN UUN=D0:GOSUB *NORMAL.PROB.HGF

1790 IF PROB$="X" THEN UUX=D0:GOSUB *CHI2.PROB.HGF

1800 IF PROB$="T" THEN UUT=D0:GOSUB *T.PROB.HGF

1810 IF PROB$="F" THEN UUF=D0:GOSUB *F.PROB.HGF

1820 IF PROB$="NX" THEN UUX=D0:GOSUB *NCHI2.PROB.HGF

1830 IF PROB$="NT" THEN UUT=D0:GOSUB *NT.PROB.HGF

1840 IF PROB$="NF" THEN UUF=D0:GOSUB *NF.PROB.HGF

1850 '

1860 F=P0-Q0:G=P1

1870 '

1880 F2=F:D2=D0

1890 IF F1=0 THEN 1920

1900　GOSUB *BISECTION

1910 'GOSUB *DAMPING

1920 '

1930 H=P2:DI=G*G-F*H*2

1940 IF DI<0 THEN DELTA=-F/G:GOTO 1970

1950 IF G<0 THEN DELTA=F*2/(-G+SQR(DI)) ELSE DELTA=F*2/(-G-SQR(DI))

1960 '

1970 DD=D0+DELTA

1980 '

1990 EPSL=.00001

2000　IF ABS(DD-D0)<EPSL THEN SW=1

2010 'IF ABS((DD-D0)/D0)<EPSL THEN SW=1

2020 F1=F2:D1=D2

2030 RETURN

2040 '

2050 *BISECTION:

2060 IF SGN(F1)=SGN(F2) THEN RETURN

2070　 DD=(D1+D2)/2:' bisection

2080　 DD=(F1*D2-F2*D1)/(F1-F2):' regula falsi

2090 RETURN 1980

2100 '

2110 *DAMPING:

2120 IF ABS(F1)>ABS(F2) THEN RETURN

2130　 DD=D0-DELTA/2:' damping

2140 RETURN 1980

2150 '

2160 ' ** NORMAL DISTRIBUTION **

2170 '

2180 *NORMAL.PROB.HGF:

2190 PI#=3.141592653589732#

2200 A0=1:C0=1.5

2210 A=A0:C=C0

2220 Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF :G1#=G#

2230 P0=.5#+G1#*UUN*EXP(-UUN*UUN/2)/SQR(PI#*2)

2240 '

2250 A=A0:C=C0+1

2260 Z=UUN*UUN/2:'GOSUB *CONFLUENT.HGF :G2#=G#*(A0/C0-1)+G1#

2270 A=A0:B=B0-1:C=C0:G0#=1#

2280 G#=(G1#-G0#)*C0/Z

2290 G2#=G#*(A0/C0-1)+G1#

2300 P10=G1#*(1-UUN*UUN)+G2#*UUN*UUN

2310 P1=P10*EXP(-UUN*UUN/2)/SQR(PI#*2)

2320 '

2330 G3#=G1#*A0/Z-G2#*(C0+Z)/Z

2340 P20=-P10*UUN+G2#*UUN*(1-UUN*UUN)-G1#*UUN*2+G3#*UUN*UUN*UUN+G2#*UUN*2

2350 P2=P20*EXP(-UUN*UUN/2)/SQR(PI#*2)

2360 RETURN

2370 '

2380 ' ** CHI-SQUARE DISTRIBUTION **

2390 '

2400 *CHI2.PROB.HGF:

2410 GOSUB *CHI2.PROB.HGF.SUB

2420 RETURN

2430 '

2440 ' ** CHI-SQUARE DISTRIBUTION **

2450 '

2460 *CHI2.PROB.HGF.SUB:

2470 A0=1:C0=1+DF/2

2480 A=A0:C=C0

2490 Z=UUX/2:GOSUB *CONFLUENT.HGF :G1#=G#

2500 Z=1+DF/2:GOSUB *LOG.GAMMA:G2#=G#

2510 Z=UUX/2

2520　P0=G1#/EXP(G2#)*Z^(DF/2)*EXP(-UUX/2)

2530 '

2540 Z=UUX/2

2550 A=A0-1:C=C0:G0#=1#

2560 G#=(G1#-G0#)*C0/Z

2570 G3#=G#*(A0/C0-1)+G1#

2580　P10=Z^(DF/2)*(DF/2/Z*G1#/2-G1#/2+G3#/2)

2590　P1=P10/EXP(G2#)*EXP(-Z)

2600 '

2610 G4#=G1#*A0/Z-G3#*(C0+Z)/Z

2620　P20=-P10/2+Z^(DF/2)*((DF/2*(DF/2-1)/Z/Z-DF/2/Z)*G1#/4+(DF/2/Z-1)*G3#/4+DF/2/Z*G3#/4+G4#/4)

2630　P2=P20/EXP(G2#)*EXP(-UUX/2)

2640 RETURN

2650 '

2660 ' ** NON-CENTRAL CHI-SQUARE DISTRIBUTION **

2670 '

2680 *NCHI2.PROB.HGF:

2690 GOSUB *CHI2.PROB.HGF.SUB

2700 GOSUB *POISSON.MIX.CHI2

2710 RETURN

2720 '

2730 ' *** ポアソン混合 ***

2740 '

2750 *POISSON.MIX.CHI2:

2760 EPSL=1E-09

2770 Z=UUX/2:P0#=P0:P1#=P1:P2#=P2:H1#=G1#

2780 P0#=EXP(-LAMBDA/2)*P0#

2790 P1#=EXP(-LAMBDA/2)*P1#

2800 P2#=EXP(-LAMBDA/2)*P2#

2810 SS#=0:S0#=P0#:S1#=P1#:S2#=P2#

2820 J=1

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

2840　 A=A0:C=C0+J-1

2850　 H2#=C/Z*(H1#-1#)

2860　 H3#=H2#*(A/C-1)+H1#

2870　 H4#=(C+1)/Z*(H2#-1#)*(A/(C+1)-1)+H2#

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

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

2900　 '

2910　 O#=H2#:I#=H1#

2920　 P0#=P0#*LAMBDA/2/J/C*O#/I#*Z

2930　 '

2940　 DFO=DF+J*2:DFI=DF+(J-1)*2

2950　 O#=DFO/2/Z*H2#/2-H2#/2+H4#/2

2960　 I#=DFI/2/Z*H1#/2-H1#/2+H3#/2

2970　 P1#=P1#*LAMBDA/2/J/C*O#/I#*Z

2980　 '

2990　 DFO=DF+J*2:DFI=DF+(J-1)*2

3000　 O#=DFO/2/Z*H2#/2-H2#/2+H4#/2

3010　 O#=-O#/2+(DFO/2*(DFO/2-1)/Z/Z-DFO/2/Z)*H2#/4+(DFO/2/Z-1)*H4#/4+DFO/2/Z*H4#/4+H6#/4

3020　 I#=DFI/2/Z*H1#/2-H1#/2+H3#/2

3030　 I#=-I#/2+(DFI/2*(DFI/2-1)/Z/Z-DFI/2/Z)*H1#/4+(DFI/2/Z-1)*H3#/4+DFI/2/Z*H3#/4+H5#/4

3040　 P2#=P2#*LAMBDA/2/J/C*O#/I#*Z

3050　 '

3060　 H1#=H2#

3070　 SS#=S0#

3080　 S0#=S0#+P0#

3090　 S1#=S1#+P1#

3100　 S2#=S2#+P2#

3110　 J=J+1

3120 WEND

3130 P0=S0#:P1=S1#:P2=S2#

3140 RETURN

3150 '

3160 ' ** T-DISTRIBUTION **

3170 '

3180 *T.PROB.HGF:

3190 DFX=DF/2:DFY=.5:GOSUB *T.PROB.HGF.SUB

3200 P0=1-P0

3210 P0=P0/2

3220 RETURN

3230 '

3240 ' ** T-DISTRIBUTION **

3250 '

3260 *T.PROB.HGF.SUB:

3270 'DFX=DF/2:DFY=.5

3280 A0=DFX+DFY:B0=1:C0=DFX+1

3290 A=A0:B=B0:C=C0

3300 Z=DF/(DF+UUT*UUT):GOSUB *GAUSS.HGF :G1#=G#

3310 P0=G1#/DFX*Z^DFX*(1-Z)^DFY

3320 Z=DFX+DFY:GOSUB *LOG.GAMMA:G2#=G#

3330 Z=DFX　　:GOSUB *LOG.GAMMA:G3#=G#

3340 Z=DFY　　:GOSUB *LOG.GAMMA:G4#=G#

3350　P0=P0*EXP(G2#-G3#-G4#)

3360 'P0=1-P0

3370 'P0=P0/2

3380 '

3390 Z=DF/(DF+UUT*UUT)

3400 A=A0+1:B=B0-1:C=C0:G0#=1#

3410 G#=-(G1#*(C0-A0-1)+G0#*(B0-C0))/(A0-B0+1)/(1-Z)

3420 G5#=(G#-G1#)*A0/Z

3430 ZZ=-DF*UUT*2/(DF+UUT*UUT)^2

3440　P10=Z^DFX*(1-Z)^DFY*((DFX/Z-DFY/(1-Z))*G1#+G5#)

3450　P1=-P10*ZZ/DFX*EXP(G2#-G3#-G4#)/2

3460 '

3470 A=A0+1:B=B0-1:C=C0:G0#=1#

3480 G6#=-G#*A0*(C0-B0*2+(B0-A0-1)*Z)/Z/(1-Z)-G0#*A0*(B0-C0)/Z/(1-Z)-G5#*C0/Z

3490 ZZZ=-DF*2*(DF-UUT*UUT*3)/(DF+UUT*UUT)^3

3500　P20=P10*ZZZ+Z^DFX*(1-Z)^DFY*(((DFX*(DFX-1)/Z/Z-DFX/Z*DFY/(1-Z)-DFX/Z*DFY/(1-Z)+DFY*(DFY-1)/(1-Z)/(1-Z))*G1#+(DFX/Z-DFY/(1-Z))*G5#+G6#)*ZZ*ZZ)

3510　P2=-P20/DFX*EXP(G2#-G3#-G4#)/2

3520 RETURN

3530 '

3540 ' ** NON-CENTRAL T-DISTRIBUTION **

3550 '

3560 *NT.PROB.HGF:

3570 PI#=3.141592653589732#

3580 UUN=-LAMBDA:GOSUB *NORMAL.PROB.HGF:T00=P0

3590 UUN= LAMBDA:'GOSUB *NORMAL.PROB.HGF

3600 UUN= LAMBDA:TTT=1-P0-.5

3610 '

3620 DFX=DF/2:DFY=.5:GOSUB *T.PROB.HGF.SUB

3630 GOSUB *POISSON.MIX.T1:T10=P0:T11=P1:T12=P2

3640 '

3650 DFX=DF/2:DFY= 1:GOSUB *T.PROB.HGF.SUB

3660 GOSUB *POISSON.MIX.T2:T20=P0:T21=P1:T22=P2

3670 '

3680 IF UUT>0 THEN P0=T00+T10+T20:P1= T11+T21:P2= T12+T22

3690 IF UUT=0 THEN P0=T00　　+T20:P1=　　 T21:P2=　　 T22

3700 IF UUT<0 THEN P0=T00-T10+T20:P1=-T11+T21:P2=-T12+T22

3710 RETURN

3720 '

3730 ' *** ポアソン混合 ***

3740 '

3750 *POISSON.MIX.T1:

3760 EPSL=1E-09

3770 Z=DF/(DF+UUT*UUT):P0#=P0:P1#=P1:P2#=P2:H1#=G1#

3780 LAMBDA2=LAMBDA*LAMBDA

3790 P0#=EXP(-LAMBDA2/2)*P0#/2

3800 P1#=EXP(-LAMBDA2/2)*P1#/2

3810 P2#=EXP(-LAMBDA2/2)*P2#/2

3820 SS#=0:S0#=P0#:S1#=P1#:S2#=P2#

3830 J=1

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

3850　 A=A0+J-1:B=B0:C=C0

3860　 H2#=-((C-A-1)*H1#+(B-C)*1#)/(A-B+1)/(1-Z)

3870　 H3#=(H2#-H1#)*A/Z

3880　 H4#=(-((C-A-2)*H2#+(B-C)*1#)/(A-B+2)/(1-Z)-H2#)*(A+1)/Z

3890　 H5#=-H2#*A*(C-B*2+(B-A-1)*Z)/Z/(1-Z)-1#*A*(B-C)/Z/(1-Z)-H3#*C/Z

3900　 H6#=((C-A-2)*H2#+(B-C)*1#)/(A-B+2)/(1-Z)*(A+1)*(C-B*2+(B-A-2)*Z)/Z/(1-Z)-1#*(A+1)*(B-C)/Z/(1-Z)-H4#*C/Z

3910　 '

3920　 O#=H2#:I#=H1#

3930　 P0#=P0#*LAMBDA2/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

3940　 '

3950　 DFYO=DFY+J:DFYI=DFY+(J-1)

3960　 O#=(DFX/Z-DFYO/(1-Z))*H2#+H4#

3970　 I#=(DFX/Z-DFYI/(1-Z))*H1#+H3#

3980　 P1#=P1#*LAMBDA2/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

3990　 '

4000　 DFYO=DFY+J:DFYI=DFY+(J-1)

4010　 O#=(DFX/Z-DFYO/(1-Z))*H2#+H4#

4020　 O#=O#*ZZZ+((DFX*(DFX-1)/Z/Z-DFX/Z*DFYO/(1-Z)-DFX/Z*DFYO/(1-Z)+DFYO*(DFYO-1)/(1-Z)/(1-Z))*H2#+(DFX/Z-DFYO/(1-Z))*H4#+H6#)*ZZ*ZZ

4030　 I#=(DFX/Z-DFYI/(1-Z))*H1#+H3#

4040　 I#=I#*ZZZ+((DFX*(DFX-1)/Z/Z-DFX/Z*DFYI/(1-Z)-DFX/Z*DFYI/(1-Z)+DFYI*(DFYI-1)/(1-Z)/(1-Z))*H1#+(DFX/Z-DFYI/(1-Z))*H3#+H5#)*ZZ*ZZ

4050　 P2#=P2#*LAMBDA2/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

4060　 '

4070　 H1#=H2#

4080　 SS#=S0#

4090　 S0#=S0#+P0#

4100　 S1#=S1#+P1#

4110　 S2#=S2#+P2#

4120　 J=J+1

4130 WEND

4140 P0=S0#:P0=.5-P0

4150 P1=S1#:P2=S2#

4160 RETURN

4170 '

4180 ' *** ポアソン混合 ***

4190 '

4200 *POISSON.MIX.T2:

4210 EPSL=1E-09

4220 Z=DF/(DF+UUT*UUT):P0#=P0:P1#=P1:P2#=P2:H1#=G1#

4230 LAMBDA2=LAMBDA*LAMBDA

4240 P0#=EXP(-LAMBDA2/2)*P0#*LAMBDA/SQR(PI#*2)

4250 P1#=EXP(-LAMBDA2/2)*P1#*LAMBDA/SQR(PI#*2)

4260 P2#=EXP(-LAMBDA2/2)*P2#*LAMBDA/SQR(PI#*2)

4270 SS#=0:S0#=P0#:S1#=P1#:S2#=P2#

4280 J=1

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

4300　 A=A0+J-1:B=B0:C=C0

4310　 H2#=-((C-A-1)*H1#+(B-C)*1#)/(A-B+1)/(1-Z)

4320　 H3#=(H2#-H1#)*A/Z

4330　 H4#=(-((C-A-2)*H2#+(B-C)*1#)/(A-B+2)/(1-Z)-H2#)*(A+1)/Z

4340　 H5#=-H2#*A*(C-B*2+(B-A-1)*Z)/Z/(1-Z)-1#*A*(B-C)/Z/(1-Z)-H3#*C/Z

4350　 H6#=((C-A-2)*H2#+(B-C)*1#)/(A-B+2)/(1-Z)*(A+1)*(C-B*2+(B-A-2)*Z)/Z/(1-Z)-1#*(A+1)*(B-C)/Z/(1-Z)-H4#*C/Z

4360　 '

4370　 O#=H2#:I#=H1#

4380　 P0#=P0#*LAMBDA2/2/(J+.5)*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

4390　 '

4400　 DFYO=DFY+J:DFYI=DFY+(J-1)

4410　 O#=(DFX/Z-DFYO/(1-Z))*H2#+H4#

4420　 I#=(DFX/Z-DFYI/(1-Z))*H1#+H3#

4430　 P1#=P1#*LAMBDA2/2/(J+.5)*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

4440　 '

4450　 DFYO=DFY+J:DFYI=DFY+(J-1)

4460　 O#=(DFX/Z-DFYO/(1-Z))*H2#+H4#

4470　 O#=O#*ZZZ+((DFX*(DFX-1)/Z/Z-DFX/Z*DFYO/(1-Z)-DFX/Z*DFYO/(1-Z)+DFYO*(DFYO-1)/(1-Z)/(1-Z))*H2#+(DFX/Z-DFYO/(1-Z))*H4#+H6#)*ZZ*ZZ

4480　 I#=(DFX/Z-DFYI/(1-Z))*H1#+H3#

4490　 I#=I#*ZZZ+((DFX*(DFX-1)/Z/Z-DFX/Z*DFYI/(1-Z)-DFX/Z*DFYI/(1-Z)+DFYI*(DFYI-1)/(1-Z)/(1-Z))*H1#+(DFX/Z-DFYI/(1-Z))*H3#+H5#)*ZZ*ZZ

4500　 P2#=P2#*LAMBDA2/2/(J+.5)*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

4510　 '

4520　 H1#=H2#

4530　 SS#=S0#

4540　 S0#=S0#+P0#

4550　 S1#=S1#+P1#

4560　 S2#=S2#+P2#

4570　 J=J+1

4580 WEND

4590 P0=S0#:P0=TTT-P0

4600 P1=S1#:P2=S2#

4610 RETURN

4620 '

4630 ' ** F-DISTRIBUTION **

4640 '

4650 *F.PROB.HGF:

4660 GOSUB *F.PROB.HGF.SUB

4670 P0=1-P0

4680 RETURN

4690 '

4700 ' ** F-DISTRIBUTION **

4710 '

4720 *F.PROB.HGF.SUB:

4730 DFX=DF2/2:DFY=DF1/2

4740 A0=DFX+DFY:B0=1:C0=DFX+1

4750 A=A0:B=B0:C=C0

4760 Z=DF2/(DF2+DF1*UUF):GOSUB *GAUSS.HGF :G1#=G#

4770 P0=G1#/DFX*Z^DFX*(1-Z)^DFY

4780 Z=DFX+DFY:GOSUB *LOG.GAMMA:G2#=G#

4790 Z=DFX　　:GOSUB *LOG.GAMMA:G3#=G#

4800 Z=DFY　　:GOSUB *LOG.GAMMA:G4#=G#

4810　P0=P0*EXP(G2#-G3#-G4#)

4820 'P0=1-P0

4830 '

4840 Z=DF2/(DF2+DF1*UUF)

4850 A=A0+1:B=B0-1:C=C0:G0#=1#

4860 G#=-(G1#*(C0-A0-1)+G0#*(B0-C0))/(A0-B0+1)/(1-Z)

4870 G5#=(G#-G1#)*A0/Z

4880 ZZ=-DF2*DF1/(DF2+DF1*UUF)^2

4890　P10=Z^DFX*(1-Z)^DFY*((DFX/Z-DFY/(1-Z))*G1#+G5#)

4900　P1=-P10*ZZ/DFX*EXP(G2#-G3#-G4#)

4910 '

4920 A=A0+1:B=B0-1:C=C0:G0#=1#

4930 G6#=-G#*A0*(C0-B0*2+(B0-A0-1)*Z)/Z/(1-Z)-G0#*A0*(B0-C0)/Z/(1-Z)-G5#*C0/Z

4940 ZZZ=-DF2*DF1*DF1*2/(DF2+DF1*UUF)^3

4950　P20=P10*ZZZ+Z^DFX*(1-Z)^DFY*(((DFX*(DFX-1)/Z/Z-DFX/Z*DFY/(1-Z)-DFX/Z*DFY/(1-Z)+DFY*(DFY-1)/(1-Z)/(1-Z))*G1#+(DFX/Z-DFY/(1-Z))*G5#+G6#)*ZZ*ZZ)

4960　P2=-P20/DFX*EXP(G2#-G3#-G4#)

4970 RETURN

4980 '

4990 ' ** NON-CENTRAL F-DISTRIBUTION **

5000 '

5010 *NF.PROB.HGF:

5020 GOSUB *F.PROB.HGF.SUB

5030 GOSUB *POISSON.MIX.F

5040 RETURN

5050 '

5060 ' *** ポアソン混合 ***

5070 '

5080 *POISSON.MIX.F:

5090 EPSL=1E-09

5100 Z=DF2/(DF2+DF1*UUF):P0#=P0:P1#=P1:P2#=P2:H1#=G1#

5110 P0#=EXP(-LAMBDA/2)*P0#

5120 P1#=EXP(-LAMBDA/2)*P1#

5130 P2#=EXP(-LAMBDA/2)*P2#

5140 SS#=0:S0#=P0#:S1#=P1#:S2#=P2#

5150 J=1

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

5170　 A=A0+J-1:B=B0:C=C0

5180　 H2#=-((C-A-1)*H1#+(B-C)*1#)/(A-B+1)/(1-Z)

5190　 H3#=(H2#-H1#)*A/Z

5200　 H4#=(-((C-A-2)*H2#+(B-C)*1#)/(A-B+2)/(1-Z)-H2#)*(A+1)/Z

5210　 H5#=-H2#*A*(C-B*2+(B-A-1)*Z)/Z/(1-Z)-1#*A*(B-C)/Z/(1-Z)-H3#*C/Z

5220　 H6#=((C-A-2)*H2#+(B-C)*1#)/(A-B+2)/(1-Z)*(A+1)*(C-B*2+(B-A-2)*Z)/Z/(1-Z)-1#*(A+1)*(B-C)/Z/(1-Z)-H4#*C/Z

5230　 '

5240　 O#=H2#:I#=H1#

5250　 P0#=P0#*LAMBDA/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

5260　 '

5270　 DFYO=DFY+J:DFYI=DFY+(J-1)

5280　 O#=(DFX/Z-DFYO/(1-Z))*H2#+H4#

5290　 I#=(DFX/Z-DFYI/(1-Z))*H1#+H3#

5300　 P1#=P1#*LAMBDA/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

5310　 '

5320　 DFYO=DFY+J:DFYI=DFY+(J-1)

5330　 O#=(DFX/Z-DFYO/(1-Z))*H2#+H4#

5340　 O#=O#*ZZZ+((DFX*(DFX-1)/Z/Z-DFX/Z*DFYO/(1-Z)-DFX/Z*DFYO/(1-Z)+DFYO*(DFYO-1)/(1-Z)/(1-Z))*H2#+(DFX/Z-DFYO/(1-Z))*H4#+H6#)*ZZ*ZZ

5350　 I#=(DFX/Z-DFYI/(1-Z))*H1#+H3#

5360　 I#=I#*ZZZ+((DFX*(DFX-1)/Z/Z-DFX/Z*DFYI/(1-Z)-DFX/Z*DFYI/(1-Z)+DFYI*(DFYI-1)/(1-Z)/(1-Z))*H1#+(DFX/Z-DFYI/(1-Z))*H3#+H5#)*ZZ*ZZ

5370　 P2#=P2#*LAMBDA/2/J*(DFX+DFY+J-1)/(DFY+J-1)*O#/I#*(1-Z)

5380　 '

5390　 H1#=H2#

5400　 SS#=S0#

5410　 S0#=S0#+P0#

5420　 S1#=S1#+P1#

5430　 S2#=S2#+P2#

5440　 J=J+1

5450 WEND

5460 P0=S0#:P0=1-P0

5470 P1=S1#:P2=S2#

5480 RETURN

5490 '

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

5510 '

5520 *CONFLUENT.HGF:' [Kummer]

5530 EPSL=1E-09

5540 G0#=0:G#=1

5550 T#=1:K=1

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

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

5580　 G0#=G#

5590　 G#=G#+T#

5600　 K=K+1

5610　 A=A+1

5620　 C=C+1

5630 WEND

5640 G=G#

5650 RETURN

5660 '

5670 ' *** ガウス型超幾何関数 ***

5680 '

5690 *GAUSS.HGF:

5700 EPSL=1E-09

5710 G0#=0:G#=1

5720 T#=1:K=1

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

5740　 T#=T#*A*B/C*Z/K

5750　 G0#=G#

5760　 G#=G#+T#

5770　 K=K+1

5780　 A=A+1

5790　 B=B+1

5800　 C=C+1

5810 WEND

5820 G=G#

5830 RETURN

5840 '

5850 ' *** 合流型超幾何関数の連分数展開 ***

5860 '

5870 *CONFLUENT.HGF1:' 1F1(1,c,x)

5880 NZ=50

5890 G#=0

5900 FOR ID=NZ TO 1 STEP -1

5910　 EVEN#=ID/(C+2*ID-1)/(C+2*ID-2)

5920　 ODD#=-(C+ID-1)/(C+2*ID)/(C+2*ID-1)

5930　 G#=ODD#*Z/(1+G#)

5940　 G#=EVEN#*Z/(1+G#)

5950 NEXT ID

5960 G#=-Z/C/(1+G#)

5970 G#=1/(1+G#)

5980 RETURN

5990 '

6000 ' *** ガウス型超幾何関数の連分数展開 ***

6010 '

6020 *GAUSS.HGF1:' 2F1(a,1,c,x)

6030 NZ=50

6040 G#=0

6050 FOR ID=NZ TO 1 STEP -1

6060　 ODD#=-(A+ID-1)*(C+ID-2)/(C+2*ID-2)/(C+2*ID-3)

6070　 EVEN#=-ID*(C-A+ID-1)/(C+2*ID-1)/(C+2*ID-2)

6080　 G#=EVEN#*Z/(1+G#)

6090　 G#=ODD#*Z/(1+G#)

6100 NEXT ID

6110 G#=1/(1+G#)

6120 RETURN

6130 '

6140 ' *** 対数ガンマ関数 (HASTINGSの多項式近似) ***

6150 '

6160 *LOG.GAMMA:

6170 A1#=-.577191652#

6180 A2#= .988205891#

6190 A3#=-.897056937#

6200 A4#= .918206857#

6210 A5#=-.756704078#

6220 A6#= .482199394#

6230 A7#=-.193527818#

6240 A8#= .035868343#

6250 X=Z:G#=0

6260 IF X<1 THEN G#=-LOG(X):GOTO 6290

6270 WHILE X>2:X=X-1:G#=G#+LOG(X):WEND

6280 X=X-1

6290 G#=G#+LOG(1+X*(A1#+X*(A2#+X*(A3#+X*(A4#+X*(A5#+X*(A6#+X*(A7#+X*A8#))))))))

6300 L.GAMMA=G#

6310 G=G#

6320 RETURN

6330 '

6340 ' ** NORMAL DISTRIBUTION **

6350 '

6360 *NORMAL.PERCENT:

6370 XXX=-LOG(4*PP*(1-PP))

6380 UUN=SQR(XXX*(2.0611786#-5.7262204#/(XXX+11.640595#)))

6390 IF PP>.5 THEN UUN=-UUN

6400 RETURN

6410 '

6420 ' ** T-DISTRIBUTION **

6430 '

6440 *T.PERCENT:

6450 GOSUB *NORMAL.PERCENT

6460 UUN=ABS(UUN)

6470 UU2=UUN*UUN

6480 AA1=(UU2+1)/DF/4

6490 AA2=((5*UU2+16)*UU2+3)/96/DF/DF

6500 AA3=(((3*UU2+19)*UU2+17)*UU2-15)/384/DF/DF/DF

6510 AA4=((((79*UU2+776)*UU2+1482)*UU2-1920)*W-945)/92460!/DF/DF/DF/DF

6520 AA5=(((((27*UU2+339)*UU2+930)*UU2-1782)*W-765)*W+17955)/368640!/DF/DF/DF/DF/DF

6530 UUT=UUN*(1+AA1+AA2+AA3+AA4+AA5)

6540 IF PP>.5 THEN UUT=-UUT

6550 RETURN

6560 '

6570 ' ** F-DISTRIBUTION **

6580 '

6590 *F.PERCENT:

6600 IF DF1=1 THEN DF=DF2:PP=PP/2:GOSUB *T.PERCENT　　　　　　　　　　　　　　　　　　　　　　　　:UUF=UUT^2:PP=PP*2:RETURN

6610 IF DF2=1 THEN DF=DF2:PP=(1-PP)/2:GOSUB *T.PERCENT　　　　　　　　　　　　　　　　　　　　　　:UUF=1/UUT^2:PP=1-PP*2:RETURN

6620 IF DF1=2 THEN UUF=DF2/2*(PP^(-2/DF2)-1):RETURN

6630 IF DF2=2 THEN UUF=2/DF1/((1-PP)^(-2/DF1)-1):RETURN

6640 AA=2/9/DF1:BB=2/9/DF2:CC=1-AA:DD=1-BB

6650 GOSUB *NORMAL.PERCENT

6660 UUF=CC*DD+UUN*SQR(CC*CC*BB+DD*DD*AA-AA*BB*UUN*UUN)

6670 UUF=(UUF/(DD*DD-BB*UUN*UUN))^3

6680 RETURN

6690 '

6700 ' ** CHI-SQUARE DISTRIBUTION **

6710 '

6720 *CHI2.PERCENT:

6730 IF DF=1 THEN PP=PP/2:GOSUB *NORMAL.PERCENT　　　　　　　　　　　　　　　　　　　　　　　　　:UUX=UUN*UUN:PP=PP*2:RETURN

6740 IF DF=2 THEN UUX=-2*LOG(PP):RETURN

6750 GOSUB *NORMAL.PERCENT

6760 UUX=2/(9*DF)

6770 UUX=1-UUX+UUN*SQR(UUX)

6780 UUX=UUX^3*DF

6790 RETURN

6800 '

6810 ' ** NORMAL DISTRIBUTION **

6820 '

6830 *NORMAL.PROB:

6840 AA1=.196854

6850 AA2=.115194

6860 AA3=.000344

6870 AA4=.019527

6880 W=ABS(UUN)

6890 P0=1+W*(AA1+W*(AA2+W*(AA3+W*AA4)))

6900 P0=P0^4

6910 P0=1-.5/P0

6920 P0=.5+(P0-.5)*SGN(UUN)

6930 PP=1-P0

6940 RETURN

6950 '

6960 ' ** T-DISTRIBUTION **

6970 '

6980 *T.PROB:

6990 W1=(1-2/(3+8*DF))

7000 W2=DF*LOG(UUT*UUT/DF+1)

7010 UUN=W1*SQR(W2)

7020 GOSUB *NORMAL.PROB

7030 PP=1-P0

7040 RETURN

7050 '

7060 ' ** F-DISTRIBUTION **

7070 '

7080 *F.PROB:

7090 FSW=0

7100 IF UUF=0 THEN P0=0:GOTO 7210

7110 IF UUF<1 THEN UUF=1/UUF:SWAP DF1,DF2:FSW=1

7120 AA1=2/(9*DF1)

7130 AA2=2/(9*DF2)

7140 W=UUF^(1/3)

7150 W1=W+AA1-W*AA2-1

7160 W2=AA2*W*W+AA1

7170 UUN=W1/SQR(W2)

7180 IF DF2<=3 THEN UUN=UUN*(1+.08*(UUN^4)/(DF2^3))

7190 GOSUB *NORMAL.PROB

7200 IF FSW=1 THEN P0=1-P0:SWAP DF1,DF2:FSW=0

7210 PP=1-P0

7220 RETURN

7230 '

7240 ' ** CHI-SQUARE DISTRIBUTION **

7250 '

7260 *CHI2.PROB:

7270 AA1=2/(9*DF)

7280 AA2=UUX/DF

7290 W=AA2^(1/3)

7300 W1=W+AA1-1

7310 UUN=W1/SQR(AA1)

7320 GOSUB *NORMAL.PROB

7330 PP=1-P0

7340 RETURN

7350 '

7360 ' ** NON-CENTRAL T-DISTRIBUTION **

7370 '

7380 *NT.PERCENT1:

7390 GOSUB *NORMAL.PERCENT1

7400 D=1-1/DF/4+1/DF/DF/32

7410 AA1=D*D

7420 AA2=D*D-UUN*UUN*(1-D*D)

7430 AA3=LAMBDA*SQR(AA1)+SGN(UUN)*SQR(LAMBDA*LAMBDA*AA1-(LAMBDA*LAMBDA-UUN*UUN)*AA2)

7440 UUT=AA3/AA2

7450 UU=UUT

7460 RETURN

7470 '

7480 '*NT.PERCENT1:

7490 GOSUB *NORMAL.PERCENT1

7500 AA1=1-1/DF/2

7510 AA2=1-1/DF/2-UUN*UUN/DF/2

7520 AA3=LAMBDA*SQR(AA1)+SGN(UUN)*SQR(LAMBDA*LAMBDA*AA1-(LAMBDA*LAMBDA-UUN*UUN)*AA2)

7530 UUT=AA3/AA2

7540 UU=UUT

7550 RETURN

7560 '

7570 ' ** NON-CENTRAL F-DISTRIBUTION **

7580 '

7590 *NF.PERCENT1:

7600 DF0=DF1

7610 DF1=DF0+LAMBDA*LAMBDA/(DF0+LAMBDA*2)

7620 C=1+LAMBDA/DF0

7630 '

7640 GOSUB *F.PERCENT1

7650 UUF=C*UUF

7660 UU=UUF:DF1=DF0

7670 RETURN

7680 '

7690 ' **　NON-CENTRAL CHI-SQUARE DISTRIBUTION **

7700 '

7710 *NCHI2.PERCENT1:

7720 DF0=DF

7730 DF=DF0+LAMBDA*LAMBDA/(DF0+LAMBDA*2)

7740 '

7750 GOSUB *NORMAL.PERCENT1

7760 AA1=2/(9*DF)

7770 W=1-AA1+UUN*SQR(AA1)

7780 UUX=(DF0+LAMBDA)*(W^3)

7790 UU=UUX:DF=DF0

7800 RETURN

7810 '

7820 '*NCHI2.PERCENT1:

7830 DF0=DF

7840 DF=DF0+LAMBDA*LAMBDA/(DF0+LAMBDA*2)

7850 C=1+LAMBDA/(DF0+LAMBDA)

7860 '

7870 GOSUB *CHI2.PERCENT1

7880 AA1=LAMBDA*LAMBDA/(DF0+LAMBDA*2)^2/(DF+2)/(DF+4)*DF/3

7890 AA2=AA1*(DF+4)*2

7900 AA3=AA1*(DF+2)*(DF+4)

7910 UUX=C*UUX*(1-UUX*UUX*AA1+UUX*AA2-AA3)

7920 UU=UUX:DF=DF0

7930 RETURN

7940 '

7950 ' ** NON-CENTRAL T-DISTRIBUTION **

7960 '

7970 *NT.PROB:

7980 D=1-1/DF/4+1/DF/DF/32

7990 AA1=LAMBDA-UUT*D

8000 AA2=SQR(1+UUT*UUT*(1-D*D))

8010 UUN=AA1/AA2

8020 GOSUB *NORMAL.PROB

8030 P0=1-P0

8040 RETURN

8050 '

8060 '*NT.PROB:

8070 AA1=LAMBDA-UUT*SQR(1-1/DF/2)

8080 AA2=SQR(1+UUT*UUT/DF/2)

8090 UUN=AA1/AA2

8100 GOSUB *NORMAL.PROB

8110 P0=1-P0

8120 RETURN

8130 '

8140 ' ** NON-CENTRAL F-DISTRIBUTION **

8150 '

8160 *NF.PROB:

8170 AA1=DF1*(UUF-1)-LAMBDA+1

8180 AA2=SQR((DF1*UUF*(1+DF1/DF2*UUF)+LAMBDA)*2)

8190 UUN=AA1/AA2

8200 GOSUB *NORMAL.PROB

8210 RETURN

8220 '

8230 ' **　NON-CENTRAL CHI-SQUARE DISTRIBUTION **

8240 '

8250 *NCHI2.PROB:

8260 AA1=UUX-LAMBDA-DF+1

8270 AA2=SQR((UUX+LAMBDA)*2)

8280 UUN=AA1/AA2

8290 GOSUB *NORMAL.PROB

8300 RETURN

8310 '

8320 ' ** NORMAL DISTRIBUTION **

8330 '

8340 *NORMAL.PERCENT1:

8350 'P0=1-PP

8360 AA1=2.30753

8370 AA2=.27061

8380 AA3=.99229

8390 AA4=.04481

8400 Q0=.5-ABS(P0-.5)

8410 W=SQR(-2*LOG(Q0))

8420 W1=AA1+AA2*W

8430 W2=1+W*(AA3+W*AA4)

8440 UUN=W-W1/W2

8450 UUN=UUN*SGN(P0-.5)

8460 UU=UUN

8470 RETURN

8480 '

8490 '*NORMAL.PERCENT1:

8500 'P0=1-PP

8510 AA1#=2.515517#

8520 AA2#=.802853#

8530 AA3#=.010328#

8540 AA4#=1.432788#

8550 AA5#=.189269#

8560 AA6#=.001308#

8570 Q0=.5-ABS(P0-.5)

8580 W=SQR(-2*LOG(Q0))

8590 W1=AA1#+W*(AA2#+W*AA3#)

8600 W2=1+W*(AA4#+W*(AA5#+W*AA6#))

8610 UUN=W-W1/W2

8620 UUN=UUN*SGN(P0-.5)

8630 UU=UUN

8640 RETURN

8650 '

8660 ' ** T-DISTRIBUTION **

8670 '

8680 *T.PERCENT1:

8690 GOSUB *NORMAL.PERCENT1

8700 W=UUN*(1+2/(1+8*DF))

8710 UUT=DF*(EXP(W*W/DF)-1)

8720 UUT=SQR(UUT)

8730 UU=UUT

8740 RETURN

8750 '

8760 ' ** F-DISTRIBUTION **

8770 '

8780 *F.PERCENT1:

8790 FSW=0

8800 IF P0<.5 THEN P0=1-P0:SWAP DF1,DF2:FSW=1

8810 GOSUB *NORMAL.PERCENT1

8820 IF DF2<=3 THEN U=UUN/(DF2^(3/4))　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　:UUN=UUN*(1.1581-.2296*U+.0042*U*U+.0027*U*U*U)

8830 AA1=2/(9*DF1)

8840 AA2=2/(9*DF2)

8850 W=UUN*UUN

8860 W1=1+AA2*(AA2-W-2)

8870 W2=AA1+AA2-AA1*AA2-1

8880 W3=1+AA1*(AA1-W-2)

8890 W4=SQR(W2*W2-W1*W3)

8900 UUF=(W4-W2)/W1

8910 UUF=UUF^3

8920 IF FSW=1 THEN UUF=1/UUF:P0=1-P0:SWAP DF1,DF2:FSW=0

8930 UU=UUF

8940 RETURN

8950 '

8960 ' ** CHI-SQUARE DISTRIBUTION **

8970 '

8980 *CHI2.PERCENT1:

8990 GOSUB *NORMAL.PERCENT1

9000 AA1=2/(9*DF)

9010 W=1-AA1+UUN*SQR(AA1)

9020 UUX=DF*(W^3)

9030 UU=UUX

9040 RETURN

　

　このプログラムの１／３は初期値ルーチンが占めています．非心分布の初期値ルーチンは簡約統計数値表（日本規格協会）と確率分布の近似（竹内啓，教育出版）を参考にして作りました．その際，本シリーズ（その１）に載せた非心分布のサブルーチンに上側確率・下側確率の混同があることが判明，その箇所を訂正して掲載しています．

　

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

　

【実行例】

　

　上のプログラムは，下側確率P0，自由度DF，非心度LAMBDAを指定することによって動かすことができます．

　

（例）

100 DF=2:LAMBDA=1

110 P0=.95:GOSUB *NCHI2.PERCENT1　 :PRINT UUX

120 P0=.95:GOSUB *NCHI2.PERCENT.HGF:PRINT UUX

130 END

　

　非心分布の計算上の問題点としては，解が真の値の近くを

　　ａ→ｂ→ａ’→ｂ’→ａ”→ｂ”→

のように振動し，一向に収束しないという状況が起こることです．この原因は，当初，無限級数Ｇ（ｘ）＝ΣｗjΣａijｘ^iＦj（ｘ）を扱っているための計算誤差が大きく，そこがボトルネックになったためだろうと思われましたが，原因はまったく別のところにありました．

　

　ｙ＝Ｇ（ｘ）は０と１の間を単調に増加する何ら変哲のない関数ですが，Ｇ’の値が０に近い，すなわちｘの大きいところでは，２次のニュートン法を用いたとしても，求める解の近傍で収束が遅延するのです．この対策としては，

　ｉ）ニュートン法の計算過程に，直線探索あるいはその簡略版であるダンピングのステップを挿入する

　ｉｉ）ニュートン法を捨て，bisection法またはregula falsi法に切り換える

ことなどが考えられました．

　

　ダンピング・ステップとは，もし新たなｘによって収束判定条件が悪化したら，Δｘ＝Δｘ／２として計算をやり直す手続きですが，１次のニュートン法ならばともかくとして，２次のニュートン法を用いているわけですから，それ程有効に働いてくれるとは思えませんし，場合によっては，収束の足を引っ張ることにもなりかねません．

　

　そこで，計算途中で

　　ａ→ｂ→ａ’→ｂ’→ａ”→ｂ”→

のように振動したならば，（ａ＋ｂ）／２としてニュートン法に戻るという単純素朴な方法で見事に収束させることができました．この方法はニュートン法とbisection法（またはregula falsi法）の折衷版と考えることができます．計算の途中まではニュートン法でせっかくうまくいっているのに，それをわざわざ全廃してまで，bisection法（またはregula falsi法）にあらためる必要はないでしょう．実際，この折衷法はニュートン法のメリットを最大限活かしながら計算しているため，ダンピング法よりもかなり有効でした．

　

（１）非心χ2分布の上側５％点の計算精度(P0=.95)

　

ν，λ　　　　 近似値　　　　　　 収束値　　　　　　　数値表

2,1　　　　　 8.5545　　　　　　 8.6422　　　　　　　8.642

2,4　　　　　 14.6591　　　　　　14.6402　　　　　　 14.641

2,16　　　　　33.3167　　　　　　33.0541　　　　　　 33.054

2,25　　　　　45.6363　　　　　　45.3121*　　　　　　45.308

4,1　　　　　 11.6691　　　　　　11.7072　　　　　　 11.707

4,4　　　　　 17.3343　　　　　　17.3093　　　　　　 17.309

4,16　　　　　35.6571　　　　　　35.427　　　　　　　35.427

4,25　　　　　47.9111　　　　　　47.6127　　　　　　 47.613

7,1　　　　　 15.9824　　　　　　16.0039　　　　　　 16.003

7,4　　　　　 22.248　　　　　　 22.2281　　　　　　 21.228

7,16　　　　　39.1613　　　　　　38.9701　　　　　　 38.970

7,25　　　　　51.3213　　　　　　51.0606　　　　　　 51.061

　

　＊の箇所は，５０回までの繰り返し計算では収束しませんでした．別の計算方法を用いると，45.3082に収束させることができました．

　

（２）非心χ2分布の下側５％点の計算精度(P0=.05)

　

ν，λ　　　　 近似値　　　　　　 収束値　　　　　　　数値表

2,1　　　　　 .170428　　　　　　.168391　　　　　　 .17

2,4　　　　　 .891717　　　　　　.645599　　　　　　 .65

2,16　　　　　6.87514　　　　　　6.32164　　　　　　 6.32

2,25　　　　　12.6707　　　　　　12.0802　　　　　　 12.08

4,1　　　　　 .910879　　　　　　.908745　　　　　　 .91

4,4　　　　　 1.9312　　　　　　 1.76501　　　　　　 1.77

4,16　　　　　8.35119　　　　　　7.88433　　　　　　 7.88

4,25　　　　　14.2598　　　　　　13.7329　　　　　　 13.73

7,1　　　　　 2.49379　　　　　　2.49371　　　　　　 2.49

7,4　　　　　 3.765　　　　　　　3.66425　　　　　　 3.66

7,16　　　　　10.6287　　　　　　10.2573　　　　　　 10.26

7,25　　　　　16.6764　　　　　　16.2267　　　　　　 16.23

　

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

　

【結語】

　

　一連のシリーズでは，シンプソンの方法など数値積分に拠らない精確な確率計算を追求してきましたが，今回（その５）では，たった１つの中心分布：Ｆ（ｘ，ν，０）から，非心分布：Ｆ（ｘ，ν，λ）の確率計算を行うコンピュータ・プログラムを紹介しました．

　

　この非心分布のプログラムは再帰的なアルゴリズムに基づいていて，キーワードは accurate, single central F, recursive です．その際，Ｆ（ｘ，ν＋２ｊ，０）を超幾何関数で表すと再帰的な関係が利用しやすくなります．

　

　これまで正規分布・χ2分布・ｔ分布・Ｆ分布の確率計算では，不完全ガンマ・不完全ベータ関数を用いる方法が行われてきましたが，超幾何関数を用いたほうが，確率計算の統一的な理解が容易になりますし，小数点自由度に対しても補間の必要がなくなるメリットは大と考えられます．

　

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