標準正規分布

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

において，累積分布関数を，

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

と定義して，この一連のコラムでは，関数列

　　Φ1(x)={1-exp(-2x^2/π)}^(1/2)

　　Φ2(x)=[1-exp(-2x^2/π){1+2(π-3)/3π^2･x^4}]^(1/2)

　　Φ3(x)=[1-exp(-2x^2/π){1+2(π-3)/3π^2･x^4-(7π^2-60π+120)/45π^3･x^6}]^(1/2)

　　Φ4(x)=[1-exp(-2x^2/π)1+2(π-3)/3π^2･x^4-(7π^2-60π+120)/45π^3･x^6+(9π^3-98π^2+420π-630)/315π^4･x^8]^(1/2)

が，Φ(x)の近似関数列になっていることをみてきた．

　

　すなわち，この順に最大絶対誤差が小さくなっていくのであるが，さらに，最大誤差を与える点がｘの大きい方へ移動することについても言及した．また，これまで触れなかった，というよりも避けて通ったのであるが，この近似関数列にはもうひとつの特徴がある．

　

　それは，誤差：Φi(x)−Φ(x)が単峰性で，山一つあるいは谷一つの曲線となるという点である．つまり，この近似関数列では

　　Φ(x)≧Φ1(x)

　　Φ(x)≦Φ2(x)

　　Φ(x)≧Φ3(x)

　　Φ(x)≦Φ4(x)

という不等式が成り立つのである．

　

　不等式が成り立つという点には捨てがたいメリットはあるのだが，最大絶対誤差を最小化させる（ミニマックス化）という観点からは，足かせになっていることも確かである．したがって，もし誤差に山と谷が現れるような近似式を得ることができたら，ミニマックス化に寄与させることができるだろうと思われる．

　

　実際，上記の不等式とは別の形式の近似式があり，Williams=山内の近似式：

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

はその例であるが，この誤差：Φ5(x)−Φ(x)は一山二谷曲線になる．

　

　そこで，今回のコラムでは，

　　Φ(x)≒{1-exp(-2x^2/π)(多項式)}^(1/2)

ではなく，別形式：

　　Φ(x)≒{1-exp(-2x^2/π)(有理式)}^(1/2)

の形の近似関数について考察することにしたい．

　

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

　

【１】パデ近似

　

　多項式近似では良い精度が得られない場合，有理式近似が良い結果を与えてくれることがある．

　

　　Φ(x)≒{1-exp(-2x^2/π)･f(x)}^(1/2)

において，多項式f(x)が

　　f(x)=Σ(-1)^k･a2k･x^(2k)

　　a2k＝1/k!2^k･4/π∫(0,π/4)(sec^2θ-π/4)^kdθ

で与えられることはわかっている．

　

　このベキ級数をあらためて，

　　f(x)=Σci･x^i

と書くことにするが，パデ近似とは，p次多項式p(x)とq次多項式q(x)からなる有理式p(x)/q(x)をとり，

　　p(x)-q(x)Σci･x^i=Σdi･x^i

の右辺の係数{di}の最初のほうが，できるだけ多く０になるとき，p(x)/q(x)をf(x)の(p,q)-パデ近似という．

　

　多項式f(x)は偶関数なので，p(x),q(x)はともに偶関数とする．ちなみに，奇関数のときはp(x)は奇関数,q(x)は偶関数とする．

　

　小生はパデ近似の経験はないのだが，その手順はプログラム化されているらしく，(p,q)=(4,2)として得られた結果が，Williams=山内の近似式：

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

であると思われる．

　

　また，ここでは述べないが，正規分布のパーセント点に関する山内の近似式もパデ近似の例なのであろう．

　

1980 '

1990 ' ** NORMAL DISTRIBUTION **

2000 '

2010 *NORMAL.PERCENT:

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

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

2040 IF PP>.5 THEN UUN=-UUN

2050 RETURN

　

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

　

【２】精度の比較

　

　以下に掲げる計算結果は，

(1)超幾何関数を使って精密計算した正規分布の下側確率の値Φ(x)

　　Φ(x)=∫(-∞,x)exp(-t^2/2)dt

　　　　　=1/2+x*EXP(-x^2/2)/√(π*2)*1F1(1,3/2;x^2/2)

と，

(2)Williamsの第一近似式Φ1(x)

(3)Williamsの第二近似式Φ2(x)

(4)Williamsの第三近似式Φ3(x)

(5)Williamsの第四近似式Φ4(x)

(6)Williams=山内の近似式Φ5(x)

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

を比較したものである．

　

x　　Ф(x)　　Ф1(x)　 Ф2(x)　 Ф3(x)　 Ф4(x)　 Ф5(x)

0.5　.691462　.691791　.691459　.691463　.691463　.691469

1.0　.841345　.843119　.841271　.841353　.841344　.841365

1.5　.933193　.936252　.932927　.93326　 .933185　.933191

2.0　.97725　 .980011　.976878　.977436　.977212　.97723

2.5　.99379　 .995301　.993534　.994025　.993718　.993793

3.0　.99865　 .999187　.998557　.998808　.998581　.998669

3.5　.999767　.999897　.99975　 .99983　 .999732　.999781

4.0　.999968　.999991　.999968　.999984　.999958　.999973

　

　ｘの値によっても異なるので一概にはいえないが，Williams=山内の近似式Φ5(x)の精度は，x=1ではФ3(x),Ф4(x)に及ばず，x=2,x=3では最も良いという結果であった．x=4ではФ2(x)の精度が最も良かったが，大ざっぱにいって，Williams=山内の近似式Φ5(x)の精度はФ4(x)以上と考えられた．

　

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

　

【３】計算プログラム

　

1000 '

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

1020 '

1030　PFILE$="scrn:"

1040 'PFILE$="b:erf.txt"

1050 OPEN PFILE$ FOR OUTPUT AS #1

1060 '

1070 FOR I=1 TO 8

1080　 UUN=I*.5

1090　 GOSUB *NORMAL.PROB.HGF:PRINT #1,P0;

1100　 GOSUB *WILLIAMS　　　 :PRINT #1,P0;

1110　 GOSUB *WILLIAMS2　　　:PRINT #1,P0;

1120　 GOSUB *WILLIAMS3　　　:PRINT #1,P0;

1130　 GOSUB *WILLIAMS4　　　:PRINT #1,P0;

1140　 GOSUB *WILLIAMS5　　　:PRINT #1,P0

1150 NEXT I

1160 CLOSE #1

1170 END

1180 '

1190 ' ** NORMAL DISTRIBUTION **

1200 '

1210 *NORMAL.PROB.HGF:

1220 PI#=3.141592653589732#

1230 A0=1:C0=1.5

1240 A=A0:C=C0

1250　Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF :G1#=G#:' [polynomial expansion]

1260 'Z=UUN*UUN/2:GOSUB *CONFLUENT.HGF1:G1#=G#:' [ continued fraction ]

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

1280 RETURN

1290 '

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

1310 '

1320 *CONFLUENT.HGF:' [Kummer]

1330 EPSL=1E-09

1340 G0#=0:G#=1

1350 T#=1:K=1

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

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

1380　 G0#=G#

1390　 G#=G#+T#

1400　 K=K+1

1410　 A=A+1

1420　 C=C+1

1430 WEND

1440 G=G#

1450 RETURN

1460 '

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

1480 '

1490 *GAUSS.HGF:

1500 EPSL=1E-09

1510 G0#=0:G#=1

1520 T#=1:K=1

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

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

1550　 G0#=G#

1560　 G#=G#+T#

1570　 K=K+1

1580　 A=A+1

1590　 B=B+1

1600　 C=C+1

1610 WEND

1620 G=G#

1630 RETURN

1640 '

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

1660 '

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

1680 NZ=20

1690 G#=0

1700 'FOR ID=NZ TO 1 STEP -1

1710 ID=NZ

1720 WHILE ID>=1

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

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

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

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

1770 ID=ID-1

1780 WEND

1790 'NEXT ID

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

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

1820 RETURN

1830 '

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

1850 '

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

1870 NZ=20

1880 G#=0

1890 'FOR ID=NZ TO 1 STEP -1

1900 ID=NZ

1910 WHILE ID>=1

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

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

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

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

1960 ID=ID-1

1970 WEND

1980 'NEXT ID

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

2000 RETURN

2010 '

2020 ' ** NORMAL DISTRIBUTION **

2030 '

2040 *WILLIAMS:

2050 UU=SQR(1-EXP(-.63662*UUN^2))/2

2060 P0=.5+UU*(SGN(UUN))

2070 RETURN

2080 '

2090 ' ** NORMAL DISTRIBUTION **

2100 '

2110 *WILLIAMS2:

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

2130 P0=.5+UU*(SGN(UUN))

2140 RETURN

2150 '

2160 ' ** NORMAL DISTRIBUTION **

2170 '

2180 *WILLIAMS3:

2190 UU=SQR(1-EXP(-.63662*UUN^2)*(1+.0095642*UUN^4-4.2405E-04*UUN^6))/2

2200 P0=.5+UU*(SGN(UUN))

2210 RETURN

2220 '

2230 ' ** NORMAL DISTRIBUTION **

2240 '

2250 *WILLIAMS4:

2260 UU=SQR(1-EXP(-.63662*UUN^2)*(1+.0095642*UUN^4-4.2405E-04*UUN^6+4.2503E-05*UUN^8))/2

2270 P0=.5+UU*(SGN(UUN))

2280 RETURN

2290 '

2300 ' ** NORMAL DISTRIBUTION **

2310 '

2320 *WILLIAMS5:

2330 UU=SQR(1-EXP(-.63662*UUN^2)*(1+UUN^4*(.0055+.0551/(UUN^2+14.4))))/2

2340 P0=.5+UU*(SGN(UUN))

2350 RETURN

　

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