（その２）で紹介した近似式：

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

は，詰まるところ，２つの独立な標準正規変数x1,x2の同時分布，すなわち，２変数正規分布の議論から，面積(2x)^2の正方形領域をそれと等しい面積の円（半径２ｘ／√π）で近似することにより得られたものである．

　

　舌足らずの点があったので補足しておくが，x1,x2が正規分布N(0,1)にしたがい，独立のとき，原点からの距離：

　　(x1^2+x2^2)^(1/2)

はレイリー分布にしたがうことになる．自由度２のχ^2分布は指数分布であるから，レイリー分布は指数分布にしたがう確率変数の平方根の分布と理解することもできよう．

　

　レイリー分布は，応用面では，２次元の標的問題（ミサイルなどの目標地点と実際の着弾地点の距離分布），いい換えればミサイルなどが目標からｒだけ離れる分布に適用されるほかに，通信工学分野（電気回路の雑音の特定の周波数について，振幅ｒと位相θとの組合せはレイリー分布に従う）など極めて重要な応用領域をもっている．

　

　また，ポアソン過程で生成された個々の点の最近接点(nearest neighbor)との距離の分布として，あるいはハザードレートを計算すると，h(x)=x/σ^2よりlinearly IFRの性質を持つ寿命分布のモデルとしても利用されている．なお，レイリー分布はχ分布であると同時にワイブル分布の１種でもある．

　

　ところで，元々の積分領域が円であるならば，たとえば正４ｎ角形においてｎ→∞とすることによっていくらでも精密に円近似することができる．しかし，ここで考えねばならないのは円ではなくて，正方形を近似することである．

　

　そのことの高次元化，すなわち，超立方体を同じ体積の超球で近似する方法を（その２）に掲げておいた．その方法は，偶数次のχ分布（多次元正規分布での原点からのユークリッド距離の確率分布で，χ^2分布の平方根分布）に帰着するわけであるが，近似式

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

は，そのｎ＝２の場合にあたっているというわけである．

　

　今回のコラムでは，高次元化した誤差関数の近似式について紹介するが，高次元化によって誤差関数の精密近似が可能になるかは甚だ疑問である．むしろ逆の結果になることが予想されるが，なぜ，このことが直観できるかという理由を述べておきたい．

　

　ｎ次元単位超立方体[-1,1]^nにおいて，単位超球が占める比率は，ｎ＝２であればπ／４(79%)であるが，ｎ＝５のときは16%に下落し，ｎ＝１０となると0.25%になる．すなわち，単位超球を超立方体中に置くと，次元が大きくなるにつれて隙間がより大きくなる．

　

　このことは，人間の直観や勘は３次元までの世界では働くものの，４次元以上の高次元についてはあてにならないという例として，しばしば取り上げられるものであるが，高次元において，超立方体内に一様分布する標本を考えるとき，低次元の場合とは対照的に大部分のデータは超球外に位置することになるのである．

　

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

　

【１】誤差関数の近似式のデザイン

　

　高次元化した近似式のデザインに必要とされる式を（その２）から拾い集めてくると，まず，不等式：

　　{Φ(x)}^n＜1/(2^(n/2-1)Γ(n/2))∫(0,r)r^(n-1)exp{-r^2/2}dr

があげられる．右辺の

　　f(x)=1/(2^(n/2-1)Γ(n/2))r^(n-1)exp{-r^2/2}

はｎ次のχ分布の密度関数，

　　F(x)=1/(2^(n/2-1)Γ(n/2))∫(0,r)r^(n-1)exp{-r^2/2}dr

はその分布関数である．

　

　ここで，

　　I[n]=∫(0,r)r^(n-1)exp{-r^2/2}dr

とおくと，部分積分より，漸化式

　　I[n]=-r^(n-2)exp{-r^2/2}+(n-2)∫(0,r)r^(n-3)exp{-r^2/2}dr

　　　　=-r^(n-2)exp{-r^2/2}+(n-2)I[n-2]　 (n>2)

が得られる．漸化式を用いると芋ヅル式に目的の値を得ることができる．

　

　また，ｎ次元超立方体を近似する超球の半径ｒは，ｎ次元単位超球の体積をｖnとすると，

　　ｖnｒ^n=(2x)^n，ｖn=π^(n/2)/Γ(n/2+1)

より

　　ｒ=2x/ｖn^(1/n)=2xπ^(-1/2){Γ(n/2+1)}^(1/n)

で与えられる．

　

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

　

【２】計算プログラム

　

　以下に掲げるプログラムは，

(1)超幾何関数を使って精密計算した正規分布の下側確率の値（1160行〜1260行）

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

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

と，

(2)Hastingsの近似式（1990行〜2130行）

(3)Cadwellの近似式≒Williams=山内の近似式（2140行〜2220行）

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

(4)Williamsの近似式，すなわち，２次元正規分布から求めた誤差関数の近似式（2260行〜2290行）

　　ΦL(x)≒1/2+1/2{1-exp(-2x^2/π)}^(1/2)

(5)４次元正規分布から求めた誤差関数の近似式（2310行〜2340行）

　　ΦL(x)≒1/2+1/2{1-exp(-2√2x^2/π)(1+2√2x^2/π)}^(1/4)

(6)６次元正規分布から求めた誤差関数の近似式（2360行〜2390行）

(7)８次元正規分布から求めた誤差関数の近似式（2410行〜2440行）

(8)１０次元正規分布から求めた誤差関数の近似式（2460行〜2490行）

の値を比較するために作成したものである．

　

　(6)(7)(8)の近似式は複雑になるので省略したが，いずれも

　　ΦL(x)≒1/2+1/2{1-exp(-2x^2/π･定数)(多項式)}^(1/n)

の形になる．すでに(5)の段階で初期モデル(4)

　　ΦL(x)≒1/2+1/2{1-exp(-2x^2/π)}^(1/2)　　　（ｘ≧０）

のエレガントさは失われてしまった感があるが，近似式の目的からすれば，やはり(4)が優れているといえるだろう．

　

　簡単なプログラムなので特別な解説は必要なかろうと思うが，2520行〜2600行で，漸化式

　　I[n]=-r^(n-2)exp{-r^2/2}+(n-2)I[n-2]　 (n>2)

の計算を行っている．

　

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 *HASTINGS　　　 :PRINT #1,P0;

1110　 GOSUB *CADWELL　　　　:PRINT #1,P0;

1120　 GOSUB *CHI

1130 NEXT I

1140 CLOSE #1

1150 END

1160 '

1170 ' ** NORMAL DISTRIBUTION **

1180 '

1190 *NORMAL.PROB.HGF:

1200 PI#=3.141592653589732#

1210 A0=1:C0=1.5

1220 A=A0:C=C0

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

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

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

1260 RETURN

1270 '

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

1290 '

1300 *CONFLUENT.HGF:' [Kummer]

1310 EPSL=1E-09

1320 G0#=0:G#=1

1330 T#=1:K=1

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

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

1360　 G0#=G#

1370　 G#=G#+T#

1380　 K=K+1

1390　 A=A+1

1400　 C=C+1

1410 WEND

1420 G=G#

1430 RETURN

1440 '

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

1460 '

1470 *GAUSS.HGF:

1480 EPSL=1E-09

1490 G0#=0:G#=1

1500 T#=1:K=1

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

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

1530　 G0#=G#

1540　 G#=G#+T#

1550　 K=K+1

1560　 A=A+1

1570　 B=B+1

1580　 C=C+1

1590 WEND

1600 G=G#

1610 RETURN

1620 '

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

1640 '

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

1660 NZ=20

1670 G#=0

1680 'FOR ID=NZ TO 1 STEP -1

1690 ID=NZ

1700 WHILE ID>=1

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

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

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

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

1750 ID=ID-1

1760 WEND

1770 'NEXT ID

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

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

1800 RETURN

1810 '

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

1830 '

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

1850 NZ=20

1860 G#=0

1870 'FOR ID=NZ TO 1 STEP -1

1880 ID=NZ

1890 WHILE ID>=1

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

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

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

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

1940 ID=ID-1

1950 WEND

1960 'NEXT ID

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

1980 RETURN

1990 '

2000 ' ** NORMAL DISTRIBUTION **

2010 '

2020 *HASTINGS:

2030 AA1=.196854

2040 AA2=.115194

2050 AA3=.000344

2060 AA4=.019527

2070 W=ABS(UUN)

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

2090 P0=P0^4

2100 P0=1-.5/P0

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

2120 PP=1-P0

2130 RETURN

2140 '

2150 ' ** NORMAL DISTRIBUTION **

2160 '

2170 *CADWELL:

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

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

2200 RETURN

2210 '

2220 ' ** NORMAL DISTRIBUTION **

2230 '

2240 *CHI:

2250 PI=3.14159

2260 R0=UUN*UUN*4/PI:C0=1

2270 R=R0^(1/2):GOSUB *CHI.SUB

2280 UU=(C0*I2)^(1/2)

2290 P0=.5+UU*(SGN(UUN))/2:PRINT #1,P0;

2300 '

2310 R0=R0*UUN*UUN*4/PI*2:C0=C0/2

2320 R=R0^(1/4):GOSUB *CHI.SUB

2330 UU=(C0*I4)^(1/4)

2340 P0=.5+UU*(SGN(UUN))/2:PRINT #1,P0;

2350 '

2360 R0=R0*UUN*UUN*4/PI*3:C0=C0/4

2370 R=R0^(1/6):GOSUB *CHI.SUB

2380 UU=(C0*I6)^(1/6)

2390 P0=.5+UU*(SGN(UUN))/2:PRINT #1,P0;

2400 '

2410 R0=R0*UUN*UUN*4/PI*4:C0=C0/6

2420 R=R0^(1/8):GOSUB *CHI.SUB

2430 UU=(C0*I8)^(1/8)

2440 P0=.5+UU*(SGN(UUN))/2:PRINT #1,P0;

2450 '

2460 R0=R0*UUN*UUN*4/PI*5:C0=C0/8

2470 R=R0^(1/10):GOSUB *CHI.SUB

2480 UU=(C0*I10)^(1/10)

2490 P0=.5+UU*(SGN(UUN))/2:PRINT #1,P0

2500 RETURN

2510 '

2520 *CHI.SUB:

2530 R2=R*R:R4=R2*R2:R6=R4*R2:R8=R6*R2

2540 ERF=EXP(-R2/2)

2550 I2=1-ERF

2560 I4=-R2*ERF+I2*2

2570 I6=-R4*ERF+I4*4

2580 I8=-R6*ERF+I6*6

2590 I10=-R8*ERF+I8*8

2600 RETURN

　

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

　

【３】計算結果

　

x　　　 (1)　　　　 (2)　　　　 (3)　　　　 (4)

0.5　　 .691462　　 .691695　　 .691459　　 .691791

1.0　　 .841345　　 .841124　　 .841271　　 .843119

1.5　　 .933193　　 .93327　　　.932927　　 .936252

2.0　　 .97725　　　.977437　　 .976878　　 .980011

2.5　　 .99379　　　.993662　　 .993534　　 .995301

3.0　　 .99865　　　.998421　　 .998557　　 .999187

3.5　　 .999767　　 .999626　　 .99975　　　.999897

4.0　　 .999968　　 .999911　　 .999968　　 .999991

　

x　　　 (5)　　　　 (6)　　　　 (7)　　　　 (8)

0.5　　 .692195　　 .69244　　　.692606　　 .692733

1.0　　 .845366　　 .846771　　 .847754　　 .848488

1.5　　 .940223　　 .942779　　 .944607　　 .946　

2.0　　 .983506　　 .985703　　 .987247　　 .988409

2.5　　 .996991　　 .997903　　 .998464　　 .998836

3.0　　 .999655　　 .999835　　 .999916　　 .999955

3.5　　 .999976　　 .999993　　 .999998　　 .999999

4.0　　 .999999　　 1　　　　　 1　　　　　 1　　　

　

　近似式の中での優劣は，

　　Cadwell=Williams=山内(3) > Hastings(2) > Williams (4)

の順であったが，さらに，

　　(4) > (5) > (6) > (7) > (8)

と続き，予想したとおり，次元が高くなるほど近似度は悪くなった．式自体が複雑になるし，そのうえ近似度が悪くなるようでは近似式の目的を逸脱してしまう．近似式としての確固たる存在理由はないといえるだろう．

　

　ところで，ｎ次元ユークリッド空間において，１辺の長さが１の立方体をｎ次元単位立方体という．その体積は１であるが，もっとも離れた２頂点を結ぶ対角線の長さはｎ次元ユークリッド空間の距離の定義から

　　√（１^2＋１^2＋・・・＋１^2）＝√ｎ

となる．

　

　したがって，次元ｎが大きくなると対角線の長さはどんどん大きくなり，ついには地球でさえ含むことができるようになる．それに対して，ｎ次元単位球はどんなに次元が高くても，長さが２より大きな線分を含むことはできない．

　

　上の計算で求めた超球（半径ｒ）と超立方体の隙間の最大距離を求めてみよう．超立方体の各面の中心までの距離はｒ−１，超立方体の頂点までの距離は√ｎ−ｒとなるので，これを表示すると，

　

n　　　　 r-1　　　　　 √n-r

2　　　　.12838　　　　.285834

4　　　　.341877　　　 .658123

6　　　　.521063　　　 .928427

8　　　　.678733　　　 1.14969

10　　　　.821266　　　 1.34101

　

　このことから，高次元になるほど超立方体とその近似超球の隙間が内外ともに大きくなっている様子を窺い知ることができる．このことが次元が高くなるほど近似度が悪くなる原因の一端となっているのである．

　

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