前々回のコラムでは，

　　１）確率分布の計算において，正規分布・χ2分布・ｔ分布・Ｆ分布の確率計算をマクローリン展開としてとらえると，収束が遅く，計算式としてはあまり実用的ではないこと，

　　２）しかし，正規分布・χ2分布・ｔ分布・Ｆ分布の分布関数は超幾何関数としてとらえることができること，

　　３）超幾何関数としてとらえなおすことによって，超幾何関数のもつ利点を活かすことができ，確率計算が見通しよくなること

を述べました．分布関数はいつも明示的な式で表されるとは限りませんから，必ずしも一筋縄では処理できないところがありますが，正規分布・χ2分布・ｔ分布・Ｆ分布の場合は，非常にラッキーなケースと考えることができるでしょう．

　

　ところで，正規分布・χ2分布・ｔ分布・Ｆ分布の片側確率やパーセント点の求め方はパッケージごとに実にさまざまで，その計算方法には，精度を優先しているもの，速度を優先しているもの，その折衷版，など少なからぬ相違がみられます．たとえば，確率変数がある値以上または以下の値をとる確率を上側確率または下側確率，ひっくるめて片側確率といいますが，

　　ａ）解析的な式に基づいて級数を足し上げる方法を採るもの（精度優先）

　　ｂ）近似式を用いているもの（速度優先）

などがあります．

　

　前者の例としては，正規分布の密度関数

　　φ（ｘ）＝１／√２πｅｘｐ（−ｘ^2／２）

の指数部分をテイラー展開（マクローリン展開）して積分すると

　　Φ(x)=1/2+1/√2π{x-x^3/2･1!･3+x^5/2^2･2!･5-･･･+(-1)^(2k+1)x^(2k+1)/2^k･k!･(2k+1)+･･･}

が得られます．収束が遅いことを厭わなければ，片側確率はこの展開式から精度よく計算することができます．

　

　χ2分布・ｔ分布・Ｆ分布では，このような展開式の引数に自由度も加わってきます．その場合，整数自由度のときは精確に計算できるのですが，自由度が整数でなく，小数点以下の端数をもつときには，

　　ａ）前後２つの整数自由度での値からの内挿補間によって求めているもの

　　ｂ）端数を切り捨てて，整数自由度に丸める

ことにより逃げているプログラムなどいろいろな対処法がみられました．要するに，解析的な式には，整数自由度でないとき精度が悪いという欠点があるのです．

　

　また，後者の例としては，正規分布におけるWilliam=山内の近似式：

　Φ(x)=1/2+1/2[1-exp(-2x^2/π){1+x^4(0.0055+0.0551/(x^2+14.4)}]^(1/2)

　

プログラム例：

1000 '

1010 ' *** 正規分布の下側確率 ***

1020 '

1030 *NORMAL.SUB:

1040 UU=SQR(1-EXP(-.63662*Z^2)*(1+.0095642*Z^4))/2

1050 P0=.5+UU*(SGN(Z)):RETURN

1060 RETURN

　

などがありますが，実際，このような近似式はとくに小さな自由度に関して信頼できませんから，純粋な数理統計学者にとってはとても気持ちの悪いことで，近似式は受け入れられないと思う人も多いかもしれません．

　

　一方，片側確率を与えたときに対応する確率変数の値をパーセント点といいます．パーセント点の計算に関しては，ほとんどの統計計算パッケージが

　　ａ）近似式を用いた計算法

を採っていて，速度優先に偏っています．

　精度優先，すなわち，

　　ｂ）ニュートン法による繰り返し計算でより精確な値を得る方法

を採用しているものは，例外的でした．

　また，小数点つきの小さい自由度に対しては，統計数値表の精度を実現できないことを理由に，

　　ｃ）数表検索と直線補間やラグランジュ補間などで補間して求める

プログラムを用意しているものもありました．

　

　今回のコラムでは，自由度が小さいときでも，小数のときでも，精確な値を保証してくれる確率分布の計算プログラムを紹介し，このプログラムから導かれた値と真値を比較してみます．

　また，超幾何関数のベキ級数展開はマクローリン展開と同じものですから，収束が遅く，収束半径なども問題になってきます．そこで，コラムの後半では，超幾何関数のベキ級数展開を連分数展開に置き換えてみます．一般に，連分数展開は速度的にも精度的にもよい計算法と信じられている節があります．実際，ベキ級数展開と連分数展開では収束速度が違ってきますが，計算精度が本当に向上するかどうかを確かめてみます．

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

　

【１】片側確率のプログラム

　

1000 '

1010 ' ** NORMAL DISTRIBUTION **

1020 '

1030 *NORMAL.PROB.HGF:

1040 PI#=3.141592653589732#

1050 'A0=.5:C0=1.5

1060 'A=A0:C=C0

1070 'Z=-UUN*UUN/2:GOSUB *CONFLUENT.HGF:G1#=G#

1080 'P0=.5+G1#*UUN/SQR(PI#*2)

1090 A0=1:C0=1.5

1100 A=A0:C=C0

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

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

1130 '

1140 'A=A0+1:C=C0+1

1150 'Z=-UUN*UUN/2:GOSUB *CONFLUENT.HGF:G2#=G#*A0/C0

1160 'P1=(G1#-G2#*UUN*UUN)/SQR(PI#*2)

1170 A=A0+1:C=C0+1

1180 Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF:G2#=G#*A0/C0

1190 P1=(G1#-G1#*UUN*UUN+G2#*UUN*UUN)*EXP(-UUN*UUN/2)/SQR(PI#*2)

1200 RETURN

1210 '

1220 ' ** CHI-SQUARE DISTRIBUTION **

1230 '

1240 *CHI2.PROB.HGF:

1250 'A0=DF/2:C0=1+DF/2

1260 'A=A0:C=C0

1270 'Z=-UUX/2:GOSUB *CONFLUENT.HGF:G1#=G#

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

1290 'P0=G1#/G2#*(UUX/2)^(DF/2)

1300 'P0=LOG(G1#)-G2#+LOG(UUX/2)*(DF/2):P0=EXP(P0)

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

1320 A=A0:C=C0

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

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

1350 'P0=G1#/G2#*(UUX/2)^(DF/2)*EXP(-UUX/2)

1360　P0=LOG(G1#)-G2#+LOG(UUX/2)*(DF/2)-UUX/2:P0=EXP(P0)

1370 '

1380 'A=A0+1:C=C0+1

1390 'Z=-UUX/2:GOSUB *CONFLUENT.HGF:G3#=G#*A0/C0

1400 'P1=(DF/2*(UUX/2)^(DF/2-1)*G1#/2-(UUX/2)^(DF/2)*G3#/2)/G2#

1410 A=A0+1:C=C0+1

1420 Z=UUX/2:GOSUB *CONFLUENT.HGF:G3#=G#*A0/C0

1430 'P1=(DF/2*(UUX/2)^(DF/2-1)*G1#/2-(UUX/2)^(DF/2)*G1#/2+(UUX/2)^(DF/2)*G3#/2)/G2#*EXP(-UUX/2)

1440　P1=LOG(DF/2*(UUX/2)^(DF/2-1)*G1#/2-(UUX/2)^(DF/2)*G1#/2+(UUX/2)^(DF/2)*G3#/2)-G2#-UUX/2:P1=EXP(P1)

1450 RETURN

1460 '

1470 ' ** T-DISTRIBUTION **

1480 '

1490 *T.PROB.HGF:

1500 DFX=DF/2:DFY=.5

1510 'A0=DFX:B0=1-DFY:C0=DFX+1

1520 'A=A0:B=B0:C=C0

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

1540 'P0=G1#/DFX*Z^DFX

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

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

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

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

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

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

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

1620 'P0=P0*G2#/G3#/G4#

1630　P0=LOG(P0)+G2#-G3#-G4#:P0=EXP(P0)

1640 P0=P0/2

1650 P0=1-P0

1660 '

1670 A=A0+1:B=B0+1:C=C0+1

1680 Z=DF/(DF+UUT*UUT):GOSUB *GAUSS.HGF:G5#=G#*A0*B0/C0

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

1700 'P1=-(DFX*Z^(DFX-1)*ZZ*G1#+Z^DFX*G5#*ZZ)/DFX*G2#/G3#/G4#/2

1710 'P1=-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ)/DFX*G2#/G3#/G4#/2

1720　P1=LOG(-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ))-LOG(DFX)+G2#-G3#-G4#-LOG(2):P1=EXP(P1)

1730 RETURN

1740 '

1750 ' ** F-DISTRIBUTION **

1760 '

1770 *F.PROB.HGF:

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

1790 'A0=DFX:B0=1-DFY:C0=DFX+1

1800 'A=A0:B=B0:C=C0

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

1820 'P0=G1#/DFX*Z^DFX

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

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

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

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

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

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

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

1900 'P0=P0*G2#/G3#/G4#

1910　P0=LOG(P0)+G2#-G3#-G4#:P0=EXP(P0)

1920 P0=1-P0

1930 '

1940 A=A0+1:B=B0+1:C=C0+1

1950 Z=DF2/(DF2+DF1*UUF):GOSUB *GAUSS.HGF:G5#=G#*A0*B0/C0

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

1970 'P1=-(DFX*Z^(DFX-1)*ZZ*G1#+Z^DFX*G5#*ZZ)/DFX*G2#/G3#/G4#

1980 'P1=-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ)/DFX*G2#/G3#/G4#

1990　P1=LOG(-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ))-LOG(DFX)+G2#-G3#-G4#:P1=EXP(P1)

2000 RETURN

2010 '

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

2030 '

2040 *CONFLUENT.HGF:' [Kummer]

2050 EPSL=1E-09

2060 G0#=0:G#=1

2070 T#=1:K=1

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

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

2100　 G0#=G#

2110　 G#=G#+T#

2120　 K=K+1

2130　 A=A+1

2140　 C=C+1

2150 WEND

2160 G=G#

2170 RETURN

2180 '

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

2200 '

2210 *GAUSS.HGF:

2220 EPSL=1E-09

2230 G0#=0:G#=1

2240 T#=1:K=1

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

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

2270　 G0#=G#

2280　 G#=G#+T#

2290　 K=K+1

2300　 A=A+1

2310　 B=B+1

2320　 C=C+1

2330 WEND

2340 G=G#

2350 RETURN

2360 '

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

2380 '

2390 *LOG.GAMMA:

2400 A1#=-.577191652#

2410 A2#= .988205891#

2420 A3#=-.897056937#

2430 A4#= .918206857#

2440 A5#=-.756704078#

2450 A6#= .482199394#

2460 A7#=-.193527818#

2470 A8#= .035868343#

2480 X=Z:G#=0

2490 IF X<1 THEN G#=-LOG(X):GOTO 2520

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

2510 X=X-1

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

2530 L.GAMMA=G#

2540 'G#=EXP(G#)

2550 G=G#

2560 RETURN

　

　核心となる部分，すなわち，超幾何関数値の計算は２０１０行〜２３５０行で行っていますが，ここでの計算は密度関数の数値積分をシンプソンの方法に拠らずに求めていることになります．また，新旧プログラムを同時に示していますが，前々回のコラムに掲げたプログラムとの変更点は，以下の通りです．

　

　ａ）超幾何関数

　正規分布は誤差関数と，χ2分布は不完全ガンマ関数，ｔ分布やＦ分布は不完全ベータ関数と密接に関係していて，正規分布，χ2分布，ベータ分布の分布関数は超幾何関数を使って以下のように表現されます．

　　Φ（ｘ）＝１／２＋１／√２πｘ1Ｆ1(1/2,3/2,-x^2/2)

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

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

　今回のプログラムでは，桁落ちを避けるために，それぞれ

　　Φ（ｘ）＝１／２＋１／√２πｘexp(-x^2/2)1F1(1,3/2,x^2/2)

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

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

と変換しました．

　

　ｂ）誤差対策としては，変数を倍精度形式にする方法が考えられますが，実際の統計計算では，入力引数の有効数字は単精度以下ですから，出力引数の桁数がいくら多くても無用の長物であり，かえって実用的ではないということになります．また，BASICの倍精度計算は非常にスピードが遅いため，全部を倍精度計算するのではなく，精度が要求される箇所だけに適宜，倍精度計算を併用しています．

　

　ｃ）１００！などの階乗計算はオーバーフローエラーを起こしますので，対数ガンマ関数を計算することによりオーバーフロー対策をしています．

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

　

【２】片側確率の計算精度

　

　正規分布の上側確率の計算精度だけを表にしますが，パーセント点の小さいところで精度がよく，大きいところでは悪くなっています．

　

正規分布の上側確率の計算精度

　

％点　　　　　 計算値　　　　　　 数値表

.2　　　　　　.42074　　　　　　 .42074

.4　　　　　　.344578　　　　　　.34458

.6　　　　　　.274253　　　　　　.27425

.8　　　　　　.211855　　　　　　.21186

1　　　　　　 .158655　　　　　　.15866

1.2　　　　　 .11507　　　　　　 .11507

1.4　　　　　 .0807567　　　　　 .080757

1.6　　　　　 .0547993　　　　　 .054799

1.8　　　　　 .0359303　　　　　 .035930

2　　　　　　 .0227501　　　　　 .022750

2.2　　　　　 .0139034　　　　　 .0103903

2.4　　　　　 8.19755E-03　　　　.0081975

2.6　　　　　 .0046612　　　　　 .0046612

2.8　　　　　 2.55501E-03　　　　.0025551

3　　　　　　 1.34981E-03　　　　.0013499

3.2　　　　　 6.87122E-04　　　　.00068714

3.4　　　　　 3.36885E-04　　　　.00033693

3.6　　　　　 1.58966E-04　　　　.00015911

3.8　　　　　 7.23004E-05　　　　.000072348

4　　　　　　 3.15905E-05　　　　.000031671

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

　

【３】パーセント点のプログラム

　

　超幾何関数の逆関数が再び超幾何関数になれば，パーセント点の計算は非常に簡単になりますが，そうは問屋が卸しません．また，ｘ→∞のとき，超幾何関数の原始関数は∞−∞になるので，広義定積分値が求められないという難点もあります．しかし，幸いなことに，超幾何関数には微分・積分してもふたたび超幾何関数になるという特性があります．

　　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)

　

　以下に掲げるパーセント点を求めるアルゴリズムでは，この性質をうまく利用して，超幾何関数の微分を行っています．そして，近似式で求まるパーセント点を初期値として，これを１次の微分まで用いたニュートン法によって補正すると，精度の良いパーセント点が得られます．

　

2570 '

2580 ' ** NORMAL DISTRIBUTION **

2590 '

2600 *NORMAL.PERCENT.HGF:

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

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

2630 UUN=X

2640 RETURN

2650 '

2660 ' ** CHI-SQUARE DISTRIBUTION **

2670 '

2680 *CHI2.PERCENT.HGF:

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

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

2710 UUX=X

2720 RETURN

2730 '

2740 ' ** T-DISTRIBUTION **

2750 '

2760 *T.PERCENT.HGF:

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

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

2790 UUT=X

2800 RETURN

2810 '

2820 ' ** F-DISTRIBUTION **

2830 '

2840 *F.PERCENT.HGF:

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

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

2870 UUF=X

2880 RETURN

2890 '

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

2910 '

2920 *NEWTON.HGF:

2930 EPSL=1E-09

2940 IMAX=50

2950 '

2960 Q0=P0

2970 D0=UU

2980 '

2990 SW=0

3000 ICNT=0

3010 WHILE SW=0

3020　 ICNT=ICNT+1

3030　 IF ICNT>IMAX THEN SW=1

3040　 DD=D0:GOSUB 3080:D0=DD

3050 WEND

3060 X=DD

3070 RETURN

3080 '

3090 ' ** 逐次計算 **

3100 '

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

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

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

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

3150 F=P0-Q0:G=P1

3160 DD=D0-F/G

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

3180 RETURN

3190 '

3200 ' ** NORMAL DISTRIBUTION **

3210 '

3220 *NORMAL.PERCENT:

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

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

3250 IF PP>.5 THEN UUN=-UUN

3260 RETURN

3270 '

3280 ' ** T-DISTRIBUTION **

3290 '

3300 *T.PERCENT:

3310 GOSUB *NORMAL.PERCENT

3320 UUN=ABS(UUN)

3330 UU2=UUN*UUN

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

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

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

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

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

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

3400 IF PP>.5 THEN UUT=-UUT

3410 RETURN

3420 '

3430 ' ** F-DISTRIBUTION **

3440 '

3450 *F.PERCENT:

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

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

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

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

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

3510 GOSUB *NORMAL.PERCENT

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

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

3540 RETURN

3550 '

3560 ' ** CHI-SQUARE DISTRIBUTION **

3570 '

3580 *CHI2.PERCENT:

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

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

3610 GOSUB *NORMAL.PERCENT

3620 UUN=ABS(UUN)

3630 UUX=2/(9*DF)

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

3650 UUX=UUX^3*DF

3660 RETURN

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

　

【４】パーセント点の計算精度

　

　最終収束値は，単精度の範囲で数値表の値とよく一致し，非常に精確な値が得られました．小数付き自由度を含む任意の自由度についての片側確率とパーセント点が簡便に求められるうえに，精度も高く，その実用上の意味は大きいと考えられます．

　

ａ）正規分布の％点の計算精度

　

上側確率　　　 近似値　　　　　　 計算値　　　　　　　数値表

.4　　　　　　.25324　　　　　　 .253347　　　　　　 .25335

.3　　　　　　.524282　　　　　　.5244　　　　　　　 .52440

.2　　　　　　.841693　　　　　　.841621　　　　　　 .84162

.1　　　　　　1.2821　　　　　　 1.28155　　　　　　 1.28155

.05　　　　　 1.64564　　　　　　1.64485　　　　　　 1.64485

.025　　　　　1.96059　　　　　　1.95996　　　　　　 1.95996

.02　　　　　 2.05426　　　　　　2.05375　　　　　　 2.05375

.01　　　　　 2.32636　　　　　　2.32635　　　　　　 2.32635

.005　　　　　2.57526　　　　　　2.57583　　　　　　 2.57583

.002　　　　　2.87696　　　　　　2.87816　　　　　　 2.87816

.001　　　　　3.08874　　　　　　3.09024　　　　　　 3.09023

　

ｂ）χ2分布の上側５％点の計算精度

　

自由度　　　　 近似値　　　　　　 計算値　　　　　　　数値表

1　　　　　　 3.84393　　　　　　3.84146　　　　　　 3.84146

2　　　　　　 5.99146　　　　　　5.99146　　　　　　 5.99146

3　　　　　　 7.77862　　　　　　7.81472　　　　　　 7.81473

4　　　　　　 9.45998　　　　　　9.48772　　　　　　 9.48773

5　　　　　　 11.0481　　　　　　11.0705　　　　　　 11.0705

6　　　　　　 12.5731　　　　　　12.5916　　　　　　 12.5916

8　　　　　　 15.4939　　　　　　15.5073　　　　　　 15.5073

10　　　　　　18.297　　　　　　 18.307　　　　　　　18.3070

12　　　　　　21.0185　　　　　　21.026　　　　　　　21.0261

15　　　　　　24.9908　　　　　　24.9958　　　　　　 24.9958

20　　　　　　31.4084　　　　　　31.4105　　　　　　 31.4104

30　　　　　　43.7744　　　　　　43.773　　　　　　　43.7730

　

ｃ）χ2分布の上側５％点の計算精度（小数自由度）

　

自由度　　　　 近似値　　　　　　 計算値　　　　　　　数値表

.5　　　　　　2.25689　　　　　　2.42024　　　　　　 2.420

.6　　　　　　2.60387　　　　　　2.7447　　　　　　　2.745

.7　　　　　　2.91994　　　　　　3.04391　　　　　　 3.044

.8　　　　　　3.21301　　　　　　3.32391　　　　　　 3.324

.9　　　　　　3.4883　　　　　　 3.5888　　　　　　　3.589

1　　　　　　 3.84393　　　　　　3.84146　　　　　　 3.841

1.1　　　　　 3.99904　　　　　　4.08401　　　　　　 4.084

1.2　　　　　 4.23904　　　　　　4.31803　　　　　　 4.318

1.3　　　　　 4.47092　　　　　　4.54478　　　　　　 4.545

1.4　　　　　 4.69583　　　　　　4.76523　　　　　　 4.765

1.5　　　　　 4.9147　　　　　　 4.98019　　　　　　 4.980

1.6　　　　　 5.12825　　　　　　5.19029　　　　　　 5.190

1.7　　　　　 5.3371　　　　　　 5.39605　　　　　　 5.396

1.8　　　　　 5.54176　　　　　　5.59793　　　　　　 5.598

1.9　　　　　 5.74264　　　　　　5.79628　　　　　　 5.796

2　　　　　　 5.99146　　　　　　5.99146　　　　　　 5.991

2.2　　　　　 6.32599　　　　　　6.37332　　　　　　 6.373

2.4　　　　　 6.70141　　　　　　6.74534　　　　　　 6.745

2.6　　　　　 7.06789　　　　　　7.10886　　　　　　 7.109

2.8　　　　　 7.42664　　　　　　7.46503　　　　　　 7.465

3　　　　　　 7.77862　　　　　　7.81472　　　　　　 7.815

3.2　　　　　 8.1246　　　　　　 8.1587　　　　　　　8.159

3.4　　　　　 8.46523　　　　　　8.49749　　　　　　 8.497

3.6　　　　　 8.80104　　　　　　8.83167　　　　　　 8.832

3.8　　　　　 9.1325　　　　　　 9.16163　　　　　　 9.162

4　　　　　　 9.45998　　　　　　9.48772　　　　　　 9.488

4.2　　　　　 9.78382　　　　　　9.81031　　　　　　 9.810

4.4　　　　　 10.1043　　　　　　10.1297　　　　　　 10.13

4.6　　　　　 10.4217　　　　　　10.446　　　　　　　10.45

4.8　　　　　 10.73638　　　　　 10.7595　　　　　　 10.76

5　　　　　　 11.0481　　　　　　11.0705　　　　　　 11.07

5.5　　　　　 11.8173　　　　　　11.8376　　　　　　 11.84

6　　　　　　 12.5731　　　　　　12.5916　　　　　　 12.59

6.5　　　　　 13.3173　　　　　　13.3343　　　　　　 13.33

7　　　　　　 14.0515　　　　　　14.0671　　　　　　 14.07

7.5　　　　　 14.7767　　　　　　14.7911　　　　　　 14.79

8　　　　　　 15.4939　　　　　　15.5073　　　　　　 15.51

8.5　　　　　 16.2039　　　　　　16.2164　　　　　　 16.22

9　　　　　　 16.9074　　　　　　16.919　　　　　　　16.92

9.5　　　　　 17.605　　　　　　 17.6157　　　　　　 17.62

10　　　　　　18.297　　　　　　 18.307　　　　　　　18.31

　

ｄ）ｔ分布の上側５％点の計算精度

　

自由度　　　　 近似値　　　　　　 計算値　　　　　　　数値表

1　　　　　　 5.64282　　　　　　6.31375　　　　　　 6.314

2　　　　　　 2.88875　　　　　　2.91999　　　　　　 2.920

3　　　　　　 2.34886　　　　　　2.35336　　　　　　 2.353

4　　　　　　 2.13136　　　　　　2.13185　　　　　　 2.132

5　　　　　　 2.01554　　　　　　2.01505　　　　　　 2.015

6　　　　　　 1.94398　　　　　　1.94318　　　　　　 1.943

7　　　　　　 1.89548　　　　　　1.89458　　　　　　 1.895

8　　　　　　 1.86048　　　　　　1.85955　　　　　　 1.860

9　　　　　　 1.83406　　　　　　1.83311　　　　　　 1.833

10　　　　　　1.8134　　　　　　 1.81246　　　　　　 1.812

12　　　　　　1.78322　　　　　　1.78229　　　　　　 1.782

15　　　　　　1.75396　　　　　　1.75305　　　　　　 1.753

20　　　　　　1.72559　　　　　　1.72472　　　　　　 1.725

30　　　　　　1.69811　　　　　　1.69726　　　　　　 1.697

50　　　　　　1.67672　　　　　　1.67591　　　　　　 1.676

100　　　　　 1.66103　　　　　　1.66023　　　　　　

　

　近似式では，自由度νの小さいところでｔ（ν，α）の誤差が大きく，信頼性に欠けます．

　

ｅ）ｔ分布の上側５％点の計算精度（小数自由度）

　

自由度　　　　 近似値　　　　　　 計算値　　　　　　　数値表

.5　　　　　　20.5678　　　　　　41.136　　　　　　　41.136

.6　　　　　　13.604　　　　　　 21.3618　　　　　　 21.362

.7　　　　　　10.0087　　　　　　13.5871　　　　　　 13.587

.8　　　　　　7.90406　　　　　　9.7866　　　　　　　9.787

.9　　　　　　6.55919　　　　　　7.64534　　　　　　 7.645

1　　　　　　 5.64282　　　　　　6.31375　　　　　　 6.314

1.1　　　　　 4.98707　　　　　　5.42384　　　　　　 5.424

1.2　　　　　 4.4994　　　　　　 4.7958　　　　　　　4.796

1.3　　　　　 4.12529　　　　　　4.33334　　　　　　 4.333

1.4　　　　　 3.8309　　　　　　 3.98105　　　　　　 3.981

1.5　　　　　 3.59426　　　　　　3.70518　　　　　　 3.705

1.6　　　　　 3.40056　　　　　　3.48415　　　　　　 3.484

1.7　　　　　 3.23956　　　　　　3.30361　　　　　　 3.304

1.8　　　　　 3.10392　　　　　　3.15373　　　　　　 3.154

1.9　　　　　 2.98831　　　　　　3.02754　　　　　　 3.028

2　　　　　　 2.88875　　　　　　2.91999　　　　　　 2.920

2.2　　　　　 2.72642　　　　　　2.74677　　　　　　 2.747

2.4　　　　　 2.60008　　　　　　2.61372　　　　　　 2.614

2.6　　　　　 2.49924　　　　　　2.50856　　　　　　 2.509

2.8　　　　　 2.41704　　　　　　2.4235　　　　　　　2.424

3　　　　　　 2.34886　　　　　　2.35336　　　　　　 2.353

3.2　　　　　 2.29147　　　　　　2.29459　　　　　　 2.295

3.4　　　　　 2.24253　　　　　　2.24467　　　　　　 2.245

3.6　　　　　 2.20034　　　　　　2.20176　　　　　　 2.202

3.8　　　　　 2.16361　　　　　　2.1645　　　　　　　2.164

4　　　　　　 2.13136　　　　　　2.13185　　　　　　 2.132

4.2　　　　　 2.10283　　　　　　2.10301　　　　　　 2.103

4.4　　　　　 2.07742　　　　　　2.07736　　　　　　 2.077

4.6　　　　　 2.05464　　　　　　2.05441　　　　　　 2.054

4.8　　　　　 2.03412　　　　　　2.03374　　　　　　 2.034

5　　　　　　 2.01554　　　　　　2.01505　　　　　　 2.015

5.5　　　　　 1.97596　　　　　　1.97528　　　　　　 1.975

6　　　　　　 1.94398　　　　　　1.94318　　　　　　 1.943

6.5　　　　　 1.9176　　　　　　 1.91674　　　　　　 1.917

7　　　　　　 1.89548　　　　　　1.89458　　　　　　 1.895

7.5　　　　　 1.87667　　　　　　1.87575　　　　　　 1.876

8　　　　　　 1.86048　　　　　　1.85955　　　　　　 1.860

8.5　　　　　 1.84641　　　　　　1.84547　　　　　　 1.845

9　　　　　　 1.83406　　　　　　1.83311　　　　　　 1.833

9.5　　　　　 1.82313　　　　　　1.82219　　　　　　 1.822

10　　　　　　1.8134　　　　　　 1.81246　　　　　　 1.812

　

ｆ）Ｆ分布の上側５％点の計算精度

　

自由度　　　　 近似値　　　　　　 計算値　　　　　　　数値表

1,1　　　　　 95.978　　　　　　 161.447　　　　　　 161.448

1,2　　　　　 17.4393　　　　　　18.5128　　　　　　 18.513

1,3　　　　　 9.99192　　　　　　10.128　　　　　　　10.128

1,4　　　　　 7.67739　　　　　　7.70865　　　　　　 7.709

1,6　　　　　 5.98612　　　　　　5.98738　　　　　　 5.987

1,8　　　　　 5.32008　　　　　　5.31766　　　　　　 5.318

1,10　　　　　4.96772　　　　　　4.9646　　　　　　　4.965

1,15　　　　　4.54626　　　　　　4.54307　　　　　　 4.543

1,30　　　　　4.17373　　　　　　4.17088　　　　　　 4.171

1,60　　　　　4.00384　　　　　　4.0012　　　　　　　4.001

　

2,1　　　　　 154.775　　　　　　199.5　　　　　　　 199.500

2,2　　　　　 19　　　　　　　　 19　　　　　　　　　19.000

2,3　　　　　 9.55209　　　　　　9.55209　　　　　　 9.552

2,4　　　　　 6.94427　　　　　　6.94427　　　　　　 6.944

2,6　　　　　 5.14325　　　　　　5.14325　　　　　　 5.143

2,8　　　　　 4.45897　　　　　　4.45897　　　　　　 4.459

2,10　　　　　4.10282　　　　　　4.10282　　　　　　 4.103

2,15　　　　　3.68232　　　　　　3.68232　　　　　　 3.682

2,30　　　　　3.31583　　　　　　3.31583　　　　　　 3.316

2,60　　　　　3.15041　　　　　　3.15041　　　　　　 3.150

　

3,1　　　　　 154.775　　　　　　215.707　　　　　　 215.707

3,2　　　　　 19.1643　　　　　　19.1643　　　　　　 19.164

3,3　　　　　 9.86395　　　　　　9.27663　　　　　　 9.277

3,4　　　　　 6.75146　　　　　　6.59138　　　　　　 6.591

3,6　　　　　 4.78438　　　　　　4.75706　　　　　　 4.757

3,8　　　　　 4.07056　　　　　　4.06618　　　　　　 4.066

3,10　　　　　3.70538　　　　　　3.70826　　　　　　 3.708

3,15　　　　　3.27918　　　　　　3.28738　　　　　　 3.287

3,30　　　　　2.91155　　　　　　2.92228　　　　　　 2.922

3,60　　　　　2.74662　　　　　　2.75807　　　　　　 2.758

　

4,1　　　　　 154.775　　　　　　224.583　　　　　　 224.583

4,2　　　　　 19.2468　　　　　　19.2468　　　　　　 19.247

4,3　　　　　 9.71849　　　　　　9.11718　　　　　　 9.117

4,4　　　　　 6.55433　　　　　　6.38823　　　　　　 6.388

4,6　　　　　 4.56561　　　　　　4.53368　　　　　　 4.534

4,8　　　　　 3.84687　　　　　　3.83785　　　　　　 3.838

4,10　　　　　3.47988　　　　　　3.47805　　　　　　 3.478

4,15　　　　　3.05221　　　　　　3.05557　　　　　　 3.056

4,30　　　　　2.68389　　　　　　2.68963　　　　　　 2.690

4,60　　　　　2.5188　　　　　　 2.52522　　　　　　 2.525

　

6,1　　　　　 154.775　　　　　　233.986　　　　　　 233.986

6,2　　　　　 19.3295　　　　　　19.3295　　　　　　 19.330

6,3　　　　　 9.55348　　　　　　8.94064　　　　　　 8.941

6,4　　　　　 6.33337　　　　　　6.16313　　　　　　 6.163

6,6　　　　　 4.31885　　　　　　4.28387　　　　　　 4.284

6,8　　　　　 3.59275　　　　　　3.58058　　　　　　 3.581

6,10　　　　　3.22228　　　　　　3.21717　　　　　　 3.217

6,15　　　　　2.79058　　　　　　2.79047　　　　　　 2.790

6,30　　　　　2.41846　　　　　　2.42052　　　　　　 2.421

6,60　　　　　2.25141　　　　　　2.25405　　　　　　 2.254

　

8,1　　　　　 154.775　　　　　　238.883　　　　　　 238.883

8,2　　　　　 19.3711　　　　　　19.371　　　　　　　19.371

8,3　　　　　 9.46287　　　　　　8.84524　　　　　　 8.845

8,4　　　　　 6.21262　　　　　　6.04104　　　　　　 6.041

8,6　　　　　 4.18272　　　　　　4.1468　　　　　　　4.147

8,8　　　　　 3.45127　　　　　　3.4381　　　　　　　3.438

8,10　　　　　3.07783　　　　　　3.07166　　　　　　 3.072

8,15　　　　　2.64209　　　　　　2.6408　　　　　　　2.641

8,30　　　　　2.26539　　　　　　2.26616　　　　　　 2.266

8,60　　　　　2.09568　　　　　　2.09697　　　　　　 2.097

　

10,1　　　　　154.775　　　　　　241.882　　　　　　 241.882

10,2　　　　　19.3958　　　　　　19.3959　　　　　　 19.396

10,3　　　　　9.4057　　　　　　 8.78552　　　　　　 8.786

10,4　　　　　6.13651　　　　　　5.96437　　　　　　 5.964

10,6　　　　　4.09623　　　　　　4.05996　　　　　　 4.060

10,8　　　　　3.36074　　　　　　3.34716　　　　　　 3.347

10,10　　　　 2.98486　　　　　　2.97824　　　　　　 2.978

10,15　　　　 2.54553　　　　　　2.54372　　　　　　 2.544

10,30　　　　 2.16441　　　　　　2.16458　　　　　　 2.165

10,60　　　　 1.99194　　　　　　1.99259　　　　　　 1.993

　

15,1　　　　　154.775　　　　　　245.949　　　　　　 245.950

15,2　　　　　19.4292　　　　　　19.4291　　　　　　 19.429

15,3　　　　　9.32601　　　　　　8.70286　　　　　　 8.703

15,4　　　　　6.03038　　　　　　5.85781　　　　　　 5.858

15,6　　　　　3.97456　　　　　　3.93806　　　　　　 3.938

15,8　　　　　3.23229　　　　　　3.21841　　　　　　 3.218

15,10　　　　 2.85201　　　　　　2.84502　　　　　　 2.845

15,15　　　　 2.40571　　　　　　2.40345　　　　　　 2.403

15,30　　　　 2.01519　　　　　　2.0148　　　　　　　2.015

15,60　　　　 1.8364　　　　　　 1.83644　　　　　　 1.836

　

30,1　　　　　154.775　　　　　　250.095　　　　　　 250.095

30,2　　　　　19.4627　　　　　　19.4624　　　　　　 19.462

30,3　　　　　9.24216　　　　　　8.61658　　　　　　 8.617

30,4　　　　　5.91852　　　　　　5.74588　　　　　　 5.746

30,6　　　　　3.8446　　　　　　 3.80816　　　　　　 3.808

30,8　　　　　3.09332　　　　　　3.07941　　　　　　 3.079

30,10　　　　 2.70664　　　　　　2.69955　　　　　　 2.700

30,15　　　　 2.24923　　　　　　2.24679　　　　　　 2.247

30,30　　　　 1.84152　　　　　　1.84087　　　　　　 1.841

30,60　　　　 1.64942　　　　　　1.64914　　　　　　 1.649

　

　ｔ分布の場合と同様，Ｆ分布の近似式でも，ν2の小さいところ，とくに４よりも小さい自由度に関して，Ｆ（ν1，ν2，α）の誤差はかなり大きくなっています．

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

　

【５】連分数展開式によるプログラム

　

1000 '

1010 ' ** NORMAL DISTRIBUTION **

1020 '

1030 *NORMAL.PROB.HGF:

1040 PI#=3.141592653589732#

1050 A0=1:C0=1.5

1060 A=A0:C=C0

1070 Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF1:G1#=G#

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

1090 '

1100 A=A0:C=C0+1

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

1120 P1=(G1#-G1#*UUN*UUN+G2#*UUN*UUN)*EXP(-UUN*UUN/2)/SQR(PI#*2)

1130 RETURN

1140 '

1150 ' ** CHI-SQUARE DISTRIBUTION **

1160 '

1170 *CHI2.PROB.HGF:

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

1190 A=A0:C=C0

1200 Z=UUX/2:GOSUB *CONFLUENT.HGF1:G1#=G#

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

1220 'P0=G1#/G2#*(UUX/2)^(DF/2)*EXP(-UUX/2)

1230　P0=LOG(G1#)-G2#+LOG(UUX/2)*(DF/2)-UUX/2:P0=EXP(P0)

1240 '

1250 A=A0:C=C0+1

1260 Z=UUX/2:GOSUB *CONFLUENT.HGF1:G3#=G#*(A0/C0-1)+G1#

1270 'P1=(DF/2*(UUX/2)^(DF/2-1)*G1#/2-(UUX/2)^(DF/2)*G1#/2+(UUX/2)^(DF/2)*G3#/2)/G2#*EXP(-UUX/2)

1280　P1=LOG(DF/2*(UUX/2)^(DF/2-1)*G1#/2-(UUX/2)^(DF/2)*G1#/2+(UUX/2)^(DF/2)*G3#/2)-G2#-UUX/2:P1=EXP(P1)

1290 RETURN

1300 '

1310 ' ** T-DISTRIBUTION **

1320 '

1330 *T.PROB.HGF:

1340 DFX=DF/2:DFY=.5

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

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

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

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

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

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

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

1420 'P0=P0*G2#/G3#/G4#

1430　P0=LOG(P0)+G2#-G3#-G4#:P0=EXP(P0)

1440 P0=P0/2

1450 P0=1-P0

1460 '

1470 A=A0+1:B=B0:C=C0

1480 Z=DF/(DF+UUT*UUT):GOSUB *GAUSS.HGF1:G5#=(G#-G1#)*A0/Z

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

1500 'P1=-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ)/DFX*G2#/G3#/G4#/2

1510　P1=LOG(-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ))-LOG(DFX)+G2#-G3#-G4#-LOG(2):P1=EXP(P1)

1520 RETURN

1530 '

1540 ' ** F-DISTRIBUTION **

1550 '

1560 *F.PROB.HGF:

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

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

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

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

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

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

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

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

1650 'P0=P0*G2#/G3#/G4#

1660　P0=LOG(P0)+G2#-G3#-G4#:P0=EXP(P0)

1670 P0=1-P0

1680 '

1690 A=A0+1:B=B0:C=C0

1700 Z=DF2/(DF2+DF1*UUF):GOSUB *GAUSS.HGF1:G5#=(G#-G1#)*A0/Z

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

1720 'P1=-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ)/DFX*G2#/G3#/G4#

1730　P1=LOG(-(DFX*Z^(DFX-1)*ZZ*(1-Z)^DFY*G1#-Z^DFX*DFY*(1-Z)^(DFY-1)*ZZ*G1#+Z^DFX*(1-Z)^DFY*G5#*ZZ))-LOG(DFX)+G2#-G3#-G4#:P1=EXP(P1)

1740 RETURN

1750 '

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

1770 '

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

1790 NZ=50

1800 G#=0

1810 FOR ID=NZ TO 1 STEP -1

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

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

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

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

1860 NEXT ID

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

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

1890 RETURN

1900 '

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

1920 '

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

1940 NZ=50

1950 G#=0

1960 FOR ID=NZ TO 1 STEP -1

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

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

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

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

2010 NEXT ID

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

2030 RETURN

　

　１７５０行〜２３５０行の合流型超幾何関数1F1(1,c,x)，ガウス型超幾何関数2F1(a,1,c,x)の連分数展開式は，一松信「特殊関数入門」森北出版を参考にしました．ベキ級数展開では，収束の様子を探りながら繰り返しを続けるか，打ち切るかを決めていましたが，収束状況の判定時間が惜しいので，連分数展開ではあらかじめ何回繰り返せば許容誤差をみたすかどうかを調べておき，無条件にその回数だけ繰り返しています（１７９０，１９４０行）．

　

　正規分布の場合，合流型超幾何関数の連分数展開はシェントンの式として知られるものに一致しますが，このプログラムでは，もっと一般性をもたせたコーディングにしてあります．また，このプログラムでは，微分値を漸化式

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

　　2F1(a+1,b,c,x)-2F1(a,b,c,x)=bx/c*2F1(a+1,b+1,c+1)

を用いて計算しています．

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

　

【６】連分数展開式による片側確率の計算精度

　

正規分布の上側確率の計算精度

　

％点　　　　　 連分数展開　　　　 ベキ級数展開　　　 数値表

.2　　　　　　.42074　　　　　　 .42074　　　　　　 .42074　　

.4　　　　　　.344578　　　　　　.344578　　　　　　.34458　　

.6　　　　　　.274253　　　　　　.274253　　　　　　.27425　　

.8　　　　　　.211855　　　　　　.211855　　　　　　.21186　　

1　　　　　　 .158655　　　　　　.158655　　　　　　.15866　　

1.2　　　　　 .11507　　　　　　 .11507　　　　　　 .11507　　

1.4　　　　　 .0807567　　　　　 .0807567　　　　　 .080757　

1.6　　　　　 .0547993　　　　　 .0547993　　　　　 .054799　

1.8　　　　　 .0359303　　　　　 .0359303　　　　　 .035930　

2　　　　　　 .0227501　　　　　 .0227501　　　　　 .022750　

2.2　　　　　 .0139036　　　　　 .0139034　　　　　 .0103903　

2.4　　　　　 8.19767E-03　　　　8.19755E-03　　　　.0081975　

2.6　　　　　 4.66144E-03　　　　.0046612　　　　　 .0046612　

2.8　　　　　 2.55531E-03　　　　2.55501E-03　　　　.0025551　

3　　　　　　 1.35034E-03　　　　1.34981E-03　　　　.0013499　

3.2　　　　　 6.87301E-04　　　　6.87122E-04　　　　.00068714

3.4　　　　　 3.38435E-04　　　　3.36885E-04　　　　.00033693

3.6　　　　　 1.65761E-04　　　　1.58966E-04　　　　.00015911

3.8　　　　　 8.40425E-05　　　　7.23004E-05　　　　.000072348

4　　　　　　 2.44975E-05　　　　3.15905E-05　　　　.000031671

　

　連分数展開式は％点の小さいところで精度がよいのですが，大きいところでは，ベキ級数展開よりもかなり悪くなっています．変数をすべて倍精度形式にして，計算項数を多くすれば精度が向上するのではないかと考えられますが，時間がかかるうえに，誤差も小さくはなりませんでした．したがって，連分数展開の精度はベキ級数展開に及ばず，精度がよいとする従来の見解には懐疑的にならざるを得ません．

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

　

【７】連分数展開式によるパーセント点の計算精度

　

ａ）正規分布の％点の計算精度

　

上側確率　　　 連分数展開　　　　 ベキ級数展開　　　　数値表

.4　　　　　　.253347　　　　　　.253347　　　　　　 .25335

.3　　　　　　.5244　　　　　　　.5244　　　　　　　 .52440

.2　　　　　　.841621　　　　　　.841621　　　　　　 .84162

.1　　　　　　1.28155　　　　　　1.28155　　　　　　 1.28155

.05　　　　　 1.64485　　　　　　1.64485　　　　　　 1.64485

.025　　　　　1.95996　　　　　　1.95996　　　　　　 1.95996

.02　　　　　 2.05375　　　　　　2.05375　　　　　　 2.05375

.01　　　　　 2.32635　　　　　　2.32635　　　　　　 2.32635

.005　　　　　2.57586　　　　　　2.57583　　　　　　 2.57583

.002　　　　　2.87816　　　　　　2.87816　　　　　　 2.87816

.001　　　　　3.09043　　　　　　3.09024　　　　　　 3.09023

　

ｂ）χ2分布の上側５％点の計算精度

　

自由度　　　　 連分数展開　　　　 ベキ級数展開　　　　数値表

1　　　　　　 3.84146　　　　　　3.84146　　　　　　 3.84146

2　　　　　　 5.99147　　　　　　5.99146　　　　　　 5.99146

3　　　　　　 7.81474　　　　　　7.81472　　　　　　 7.81473

4　　　　　　 9.48773　　　　　　9.48772　　　　　　 9.48773

5　　　　　　 11.0705　　　　　　11.0705　　　　　　 11.0705

6　　　　　　 12.5916　　　　　　12.5916　　　　　　 12.5916

8　　　　　　 15.5073　　　　　　15.5073　　　　　　 15.5073

10　　　　　　18.307　　　　　　 18.307　　　　　　　18.3070

12　　　　　　21.026　　　　　　 21.026　　　　　　　21.0261

15　　　　　　24.9958　　　　　　24.9958　　　　　　 24.9958

20　　　　　　31.4105　　　　　　31.4105　　　　　　 31.4104

30　　　　　　43.7728　　　　　　43.773　　　　　　　43.7730

　

　この連分数展開では，自由度が大きいとき，パーセント点の小さいところでの収束が遅い．

　

ｃ）χ2分布の上側５％点の計算精度（小数自由度）

　

自由度　　　　 連分数展開　　　　 ベキ級数展開　　　　数値表

.5　　　　　　2.42024　　　　　　2.42024　　　　　　 2.420

.6　　　　　　2.7447　　　　　　 2.7447　　　　　　　2.745

.7　　　　　　3.04391　　　　　　3.04391　　　　　　 3.044

.8　　　　　　3.32391　　　　　　3.32391　　　　　　 3.324

.9　　　　　　3.5888　　　　　　 3.5888　　　　　　　3.589

1　　　　　　 3.84146　　　　　　3.84146　　　　　　 3.841

1.1　　　　　 4.084　　　　　　　4.08401　　　　　　 4.084

1.2　　　　　 4.31802　　　　　　4.31803　　　　　　 4.318

1.3　　　　　 4.54478　　　　　　4.54478　　　　　　 4.545

1.4　　　　　 4.76523　　　　　　4.76523　　　　　　 4.765

1.5　　　　　 4.98019　　　　　　4.98019　　　　　　 4.980

1.6　　　　　 5.19029　　　　　　5.19029　　　　　　 5.190

1.7　　　　　 5.39605　　　　　　5.39605　　　　　　 5.396

1.8　　　　　 5.59794　　　　　　5.59793　　　　　　 5.598

1.9　　　　　 5.79629　　　　　　5.79628　　　　　　 5.796

2　　　　　　 5.99147　　　　　　5.99146　　　　　　 5.991

2.2　　　　　 6.37329　　　　　　6.37332　　　　　　 6.373

2.4　　　　　 6.7453　　　　　　 6.74534　　　　　　 6.745

2.6　　　　　 7.10886　　　　　　7.10886　　　　　　 7.109

2.8　　　　　 7.46502　　　　　　7.46503　　　　　　 7.465

3　　　　　　 7.81474　　　　　　7.81472　　　　　　 7.815

3.2　　　　　 8.15866　　　　　　8.1587　　　　　　　8.159

3.4　　　　　 8.49746　　　　　　8.49749　　　　　　 8.497

3.6　　　　　 8.83164　　　　　　8.83167　　　　　　 8.832

3.8　　　　　 9.16162　　　　　　9.16163　　　　　　 9.162

4　　　　　　 9.48773　　　　　　9.48772　　　　　　 9.488

4.2　　　　　 9.8103　　　　　　 9.81031　　　　　　 9.810

4.4　　　　　 10.1296　　　　　　10.1297　　　　　　 10.13

4.6　　　　　 10.446　　　　　　 10.446　　　　　　　10.45

4.8　　　　　 10.7596　　　　　　10.7595　　　　　　 10.76

5　　　　　　 11.0705　　　　　　11.0705　　　　　　 11.07

5.5　　　　　 11.8376　　　　　　11.8376　　　　　　 11.84

6　　　　　　 12.5916　　　　　　12.5916　　　　　　 12.59

6.5　　　　　 13.3344　　　　　　13.3343　　　　　　 13.33

7　　　　　　 14.0672　　　　　　14.0671　　　　　　 14.07

7.5　　　　　 14.7912　　　　　　14.7911　　　　　　 14.79

8　　　　　　 15.5073　　　　　　15.5073　　　　　　 15.51

8.5　　　　　 16.2163　　　　　　16.2164　　　　　　 16.22

9　　　　　　 16.919　　　　　　 16.919　　　　　　　16.92

9.5　　　　　 17.6157　　　　　　17.6157　　　　　　 17.62

10　　　　　　18.307　　　　　　 18.307　　　　　　　18.31

　

ｄ）ｔ分布の上側５％点の計算精度

　

自由度　　　　 連分数展開　　　　 ベキ級数展開　　　　数値表

1　　　　　　 6.31375　　　　　　6.31375　　　　　　 6.314

2　　　　　　 2.91999　　　　　　2.91999　　　　　　 2.920

3　　　　　　 2.35336　　　　　　2.35336　　　　　　 2.353

4　　　　　　 2.13185　　　　　　2.13185　　　　　　 2.132

5　　　　　　 2.01505　　　　　　2.01505　　　　　　 2.015

6　　　　　　 1.94318　　　　　　1.94318　　　　　　 1.943

7　　　　　　 1.89458　　　　　　1.89458　　　　　　 1.895

8　　　　　　 1.85955　　　　　　1.85955　　　　　　 1.860

9　　　　　　 1.83311　　　　　　1.83311　　　　　　 1.833

10　　　　　　1.81246　　　　　　1.81246　　　　　　 1.812

12　　　　　　1.78229　　　　　　1.78229　　　　　　 1.782

15　　　　　　1.75305　　　　　　1.75305　　　　　　 1.753

20　　　　　　1.72472　　　　　　1.72472　　　　　　 1.725

30　　　　　　1.69726　　　　　　1.69726　　　　　　 1.697

50　　　　　　1.6759　　　　　　 1.67591　　　　　　 1.676

100　　　　　 1.66023　　　　　　1.66023　　　　　　　　　

　

ｅ）ｔ分布の上側５％点の計算精度（小数自由度）

　

自由度　　　　 連分数展開　　　　 ベキ級数展開　　　　数値表

.5　　　　　　41.136　　　　　　 41.136　　　　　　　41.136

.6　　　　　　21.3618　　　　　　21.3618　　　　　　 21.362

.7　　　　　　13.5871　　　　　　13.5871　　　　　　 13.587

.8　　　　　　9.7866　　　　　　 9.7866　　　　　　　9.787　

.9　　　　　　7.64534　　　　　　7.64534　　　　　　 7.645　

1　　　　　　 6.31375　　　　　　6.31375　　　　　　 6.314　

1.1　　　　　 5.42384　　　　　　5.42384　　　　　　 5.424　

1.2　　　　　 4.7958　　　　　　 4.7958　　　　　　　4.796　

1.3　　　　　 4.33334　　　　　　4.33334　　　　　　 4.333　

1.4　　　　　 3.98105　　　　　　3.98105　　　　　　 3.981　

1.5　　　　　 3.70518　　　　　　3.70518　　　　　　 3.705　

1.6　　　　　 3.48415　　　　　　3.48415　　　　　　 3.484　

1.7　　　　　 3.30361　　　　　　3.30361　　　　　　 3.304　

1.8　　　　　 3.15373　　　　　　3.15373　　　　　　 3.154　

1.9　　　　　 3.02754　　　　　　3.02754　　　　　　 3.028　

2　　　　　　 2.91999　　　　　　2.91999　　　　　　 2.920　

2.2　　　　　 2.74677　　　　　　2.74677　　　　　　 2.747　

2.4　　　　　 2.61372　　　　　　2.61372　　　　　　 2.614　

2.6　　　　　 2.50856　　　　　　2.50856　　　　　　 2.509　

2.8　　　　　 2.4235　　　　　　 2.4235　　　　　　　2.424　

3　　　　　　 2.35336　　　　　　2.35336　　　　　　 2.353　

3.2　　　　　 2.29459　　　　　　2.29459　　　　　　 2.295　

3.4　　　　　 2.24467　　　　　　2.24467　　　　　　 2.245　

3.6　　　　　 2.20176　　　　　　2.20176　　　　　　 2.202　

3.8　　　　　 2.1645　　　　　　 2.1645　　　　　　　2.164　

4　　　　　　 2.13185　　　　　　2.13185　　　　　　 2.132　

4.2　　　　　 2.10301　　　　　　2.10301　　　　　　 2.103　

4.4　　　　　 2.07736　　　　　　2.07736　　　　　　 2.077　

4.6　　　　　 2.05441　　　　　　2.05441　　　　　　 2.054　

4.8　　　　　 2.03374　　　　　　2.03374　　　　　　 2.034　

5　　　　　　 2.01505　　　　　　2.01505　　　　　　 2.015　

5.5　　　　　 1.97528　　　　　　1.97528　　　　　　 1.975　

6　　　　　　 1.94318　　　　　　1.94318　　　　　　 1.943　

6.5　　　　　 1.91674　　　　　　1.91674　　　　　　 1.917　

7　　　　　　 1.89458　　　　　　1.89458　　　　　　 1.895　

7.5　　　　　 1.87575　　　　　　1.87575　　　　　　 1.876　

8　　　　　　 1.85955　　　　　　1.85955　　　　　　 1.860　

8.5　　　　　 1.84547　　　　　　1.84547　　　　　　 1.845　

9　　　　　　 1.83311　　　　　　1.83311　　　　　　 1.833　

9.5　　　　　 1.82219　　　　　　1.82219　　　　　　 1.822　

10　　　　　　1.81246　　　　　　1.81246　　　　　　 1.812　

　

ｆ）Ｆ分布の上側５％点の計算精度

　

自由度　　　　 連分数展開　　　　 ベキ級数展開　　　　数値表

1,1　　　　　 161.447　　　　　　161.447　　　　　　 161.448

1,2　　　　　 18.5128　　　　　　18.5128　　　　　　 18.513

1,3　　　　　 10.128　　　　　　 10.128　　　　　　　10.128

1,4　　　　　 7.70865　　　　　　7.70865　　　　　　 7.709　

1,6　　　　　 5.98738　　　　　　5.98738　　　　　　 5.987　

1,8　　　　　 5.31766　　　　　　5.31766　　　　　　 5.318　

1,10　　　　　4.9646　　　　　　 4.9646　　　　　　　4.965　

1,15　　　　　4.54307　　　　　　4.54307　　　　　　 4.543　

1,30　　　　　4.17088　　　　　　4.17088　　　　　　 4.171　

1,60　　　　　4.00119　　　　　　4.0012　　　　　　　4.001　

2,1　　　　　 199.5　　　　　　　199.5　　　　　　　 199.500

2,2　　　　　 19　　　　　　　　 19　　　　　　　　　19.000

2,3　　　　　 9.55209　　　　　　9.55209　　　　　　 9.552　

2,4　　　　　 6.94427　　　　　　6.94427　　　　　　 6.944　

2,6　　　　　 5.14325　　　　　　5.14325　　　　　　 5.143　

2,8　　　　　 4.45897　　　　　　4.45897　　　　　　 4.459　

2,10　　　　　4.10282　　　　　　4.10282　　　　　　 4.103　

2,15　　　　　3.68232　　　　　　3.68232　　　　　　 3.682　

2,30　　　　　3.31583　　　　　　3.31583　　　　　　 3.316　

2,60　　　　　3.15041　　　　　　3.15041　　　　　　 3.150　

3,1　　　　　 215.707　　　　　　215.707　　　　　　 215.707

3,2　　　　　 19.1643　　　　　　19.1643　　　　　　 19.164

3,3　　　　　 9.27663　　　　　　9.27663　　　　　　 9.277　

3,4　　　　　 6.59138　　　　　　6.59138　　　　　　 6.591　

3,6　　　　　 4.75706　　　　　　4.75706　　　　　　 4.757　

3,8　　　　　 4.06618　　　　　　4.06618　　　　　　 4.066　

3,10　　　　　3.70826　　　　　　3.70826　　　　　　 3.708　

3,15　　　　　3.28738　　　　　　3.28738　　　　　　 3.287　

3,30　　　　　2.92228　　　　　　2.92228　　　　　　 2.922　

3,60　　　　　2.75807　　　　　　2.75807　　　　　　 2.758　

　

4,1　　　　　 224.583　　　　　　224.583　　　　　　 224.583

4,2　　　　　 19.2468　　　　　　19.2468　　　　　　 19.247

4,3　　　　　 9.11718　　　　　　9.11718　　　　　　 9.117　

4,4　　　　　 6.38823　　　　　　6.38823　　　　　　 6.388　

4,6　　　　　 4.53368　　　　　　4.53368　　　　　　 4.534　

4,8　　　　　 3.83785　　　　　　3.83785　　　　　　 3.838　

4,10　　　　　3.47805　　　　　　3.47805　　　　　　 3.478　

4,15　　　　　3.05557　　　　　　3.05557　　　　　　 3.056　

4,30　　　　　2.68963　　　　　　2.68963　　　　　　 2.690　

4,60　　　　　2.52522　　　　　　2.52522　　　　　　 2.525　

　

6,1　　　　　 233.986　　　　　　233.986　　　　　　 233.986

6,2　　　　　 19.3295　　　　　　19.3295　　　　　　 19.330

6,3　　　　　 8.94064　　　　　　8.94064　　　　　　 8.941　

6,4　　　　　 6.16313　　　　　　6.16313　　　　　　 6.163　

6,6　　　　　 4.28387　　　　　　4.28387　　　　　　 4.284　

6,8　　　　　 3.58058　　　　　　3.58058　　　　　　 3.581　

6,10　　　　　3.21717　　　　　　3.21717　　　　　　 3.217　

6,15　　　　　2.79047　　　　　　2.79047　　　　　　 2.790　

6,30　　　　　2.42052　　　　　　2.42052　　　　　　 2.421　

6,60　　　　　2.25405　　　　　　2.25405　　　　　　 2.254　

　

8,1　　　　　 238.883　　　　　　238.883　　　　　　 238.883

8,2　　　　　 19.371　　　　　　 19.371　　　　　　　19.371

8,3　　　　　 8.84524　　　　　　8.84524　　　　　　 8.845　

8,4　　　　　 6.04104　　　　　　6.04104　　　　　　 6.041　

8,6　　　　　 4.1468　　　　　　 4.1468　　　　　　　4.147　

8,8　　　　　 3.4381　　　　　　 3.4381　　　　　　　3.438　

8,10　　　　　3.07166　　　　　　3.07166　　　　　　 3.072　

8,15　　　　　2.64408　　　　　　2.6408　　　　　　　2.641　

8,30　　　　　2.26616　　　　　　2.26616　　　　　　 2.266　

8,60　　　　　2.09697　　　　　　2.09697　　　　　　 2.097　

　

10,1　　　　　241.882　　　　　　241.882　　　　　　 241.882

10,2　　　　　19.3959　　　　　　19.3959　　　　　　 19.396

10,3　　　　　8.78552　　　　　　8.78552　　　　　　 8.786　

10,4　　　　　5.96437　　　　　　5.96437　　　　　　 5.964　

10,6　　　　　4.05996　　　　　　4.05996　　　　　　 4.060　

10,8　　　　　3.34716　　　　　　3.34716　　　　　　 3.347　

10,10　　　　 2.97824　　　　　　2.97824　　　　　　 2.978　

10,15　　　　 2.54372　　　　　　2.54372　　　　　　 2.544　

10,30　　　　 2.16458　　　　　　2.16458　　　　　　 2.165　

10,60　　　　 1.99259　　　　　　1.99259　　　　　　 1.993　

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

15,1　　　　　245.949　　　　　　245.949　　　　　　 245.950

15,2　　　　　19.4291　　　　　　19.4291　　　　　　 19.429

15,3　　　　　8.70286　　　　　　8.70286　　　　　　 8.703　

15,4　　　　　5.85781　　　　　　5.85781　　　　　　 5.858　

15,6　　　　　3.93806　　　　　　3.93806　　　　　　 3.938　

15,8　　　　　3.21841　　　　　　3.21841　　　　　　 3.218　

15,10　　　　 2.84502　　　　　　2.84502　　　　　　 2.845　

15,15　　　　 2.40345　　　　　　2.40345　　　　　　 2.403　

15,30　　　　 2.0148　　　　　　 2.0148　　　　　　　2.015　

15,60　　　　 1.83644　　　　　　1.83644　　　　　　 1.836　

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

30,1　　　　　250.095　　　　　　250.095　　　　　　 250.095

30,2　　　　　19.4624　　　　　　19.4624　　　　　　 19.462

30,3　　　　　8.61658　　　　　　8.61658　　　　　　 8.617　

30,4　　　　　5.74588　　　　　　5.74588　　　　　　 5.746　

30,6　　　　　3.80816　　　　　　3.80816　　　　　　 3.808　

30,8　　　　　3.07941　　　　　　3.07941　　　　　　 3.079　

30,10　　　　 2.69955　　　　　　2.69955　　　　　　 2.700　

30,15　　　　 2.24679　　　　　　2.24679　　　　　　 2.247　

30,30　　　　 1.84087　　　　　　1.84087　　　　　　 1.841　

30,60　　　　 1.64914　　　　　　1.64914　　　　　　 1.649　

　

　一般的にいって，連分数展開は計算の手間に比して精度がよくないという期待はずれの結果となりました．連分数展開にとっては悲観的な見方になりますが，高精度計算が必要なときは，ベキ級数展開のほうがよいし，高精度計算が必要なとき以外は，近似式を用いるほうが得策と考えられます．

　ただし，近似式では自由度νの小さいところでｔ（ν，α），自由度ν2の小さいところでＦ（ν1，ν2，α）の誤差がかなり大きいことを念頭に入れておく必要があります．

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