今回のコラムでは，実際にプログラムを動かして，最近接球面距離分布の分布関数の計算結果を示してみたい．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】計算プログラム
一般的にいって，
(1)連分数展開は計算の手間に比して精度がよくない
(2)高精度計算が必要なときは，ベキ級数展開のほうがよい
のですが，ベキ級数展開ではθの小さいところで収束しなっかったため，連分数展開を用いました．１１８０行で計算結果が出力されますが，Ｔ，Ｆ，Ｇはそれぞれ球面距離θ，球面上の最近接距離分布の確率密度関数，分布関数です．
1000 '
1010 ' **** nearest neighbor on S^(n-1) ****
1020 '
1030 PI=3.14159
1040 N=3 :' [dimension]
1050 NN=100:' [data no.]
1060 '
1070 DFX=1/2 :DFY=(N-1)/2
1080 Z=DFX+DFY:GOSUB *LOG.GAMMA:G2#=G#
1090 Z=DFX :GOSUB *LOG.GAMMA:G3#=G#
1100 Z=DFY :GOSUB *LOG.GAMMA:G4#=G#
1110 '
1120 FOR T=0 TO 180 STEP 10
1130 S=1:IF T>90 THEN S=-1
1140 TT=T/180*PI
1150 IF N=2 THEN F=NN/PI*EXP(-NN/PI*TT)/(1-EXP(-NN)) :G=(1-EXP(-NN/PI*TT))/(1-EXP(-NN))
1160 IF N=3 THEN F=NN*SIN(TT)/2*EXP(-NN*(1-COS(TT))/2)/(1-EXP(-NN)) :G=(1-EXP(-NN*(1-COS(TT))/2))/(1-EXP(-NN))
1170 IF N>3 THEN GOSUB *BETA.PROB.HGF.SUB
1180 PRINT T,F,G
1190 NEXT T
1200 END
1210 '
1220 ' ** BETA-DISTRIBUTION **
1230 '
1240 *BETA.PROB.HGF.SUB:
1250 'DFX=1/2 :DFY=(N-1)/2
1260 A0=DFX+DFY:B0=1:C0=DFX+1
1270 A=A0:B=B0:C=C0
1280 'Z=COS(TT)*COS(TT) :GOSUB *GAUSS.HGF :G1#=G#:' [polynomial expansion]
1290 Z=COS(TT)*COS(TT) :GOSUB *GAUSS.HGF1:G1#=G#:' [ continued fraction ]
1300 P0=G1#/DFX*Z^DFX*(1-Z)^DFY
1310 P0=P0*EXP(G2#-G3#-G4#)
1320 P1=Z^(DFX-1)*(1-Z)^(DFY-1)
1330 P1=P1*EXP(G2#-G3#-G4#)
1340 Q0=NN*(1-P0*S)/2
1350 Q1=NN*P1*COS(TT)*SIN(TT)*S
1360 F=Q1*EXP(-Q0)/(1-EXP(-NN))
1370 G=(1-EXP(-Q0))/(1-EXP(-NN))
1380 RETURN
1390 '
1400 ' *** 合流型超幾何関数 ***
1410 '
1420 *CONFLUENT.HGF:' [Kummer]
1430 EPSL=1E-09
1440 G0#=0:G#=1
1450 T#=1:K=1
1460 WHILE ABS(G0#-G#)>EPSL
1470 T#=T#*A/C*Z/K
1480 G0#=G#
1490 G#=G#+T#
1500 K=K+1
1510 A=A+1
1520 C=C+1
1530 WEND
1540 G=G#
1550 RETURN
1560 '
1570 ' *** ガウス型超幾何関数 ***
1580 '
1590 *GAUSS.HGF:
1600 EPSL=1E-09
1610 G0#=0:G#=1
1620 T#=1:K=1
1630 WHILE ABS(G0#-G#)>EPSL
1640 T#=T#*A*B/C*Z/K
1650 G0#=G#
1660 G#=G#+T#
1670 K=K+1
1680 A=A+1
1690 B=B+1
1700 C=C+1
1710 WEND
1720 G=G#
1730 RETURN
1740 '
1750 ' *** 合流型超幾何関数の連分数展開 ***
1760 '
1770 *CONFLUENT.HGF1:' 1F1(1,c,x)
1780 NZ=20
1790 G#=0
1800 'FOR ID=NZ TO 1 STEP -1
1810 ID=NZ
1820 WHILE ID>=1
1830 EVEN#=ID/(C+2*ID-1)/(C+2*ID-2)
1840 ODD#=-(C+ID-1)/(C+2*ID)/(C+2*ID-1)
1850 G#=ODD#*Z/(1#+G#)
1860 G#=EVEN#*Z/(1#+G#)
1870 ID=ID-1
1880 WEND
1890 'NEXT ID
1900 G#=-Z/C/(1#+G#)
1910 G#=1#/(1#+G#)
1920 RETURN
1930 '
1940 ' *** ガウス型超幾何関数の連分数展開 ***
1950 '
1960 *GAUSS.HGF1:' 2F1(a,1,c,x)
1970 NZ=20
1980 G#=0
1990 'FOR ID=NZ TO 1 STEP -1
2000 ID=NZ
2010 WHILE ID>=1
2020 ODD#=-(A+ID-1)*(C+ID-2)/(C+2*ID-2)/(C+2*ID-3)
2030 EVEN#=-ID*(C-A+ID-1)/(C+2*ID-1)/(C+2*ID-2)
2040 G#=EVEN#*Z/(1#+G#)
2050 G#=ODD#*Z/(1#+G#)
2060 ID=ID-1
2070 WEND
2080 'NEXT ID
2090 G#=1#/(1#+G#)
2100 RETURN
2110 '
2120 ' *** 対数ガンマ関数 (HASTINGSの多項式近似) ***
2130 '
2140 *LOG.GAMMA:
2150 A1#=-.577191652#
2160 A2#= .988205891#
2170 A3#=-.897056937#
2180 A4#= .918206857#
2190 A5#=-.756704078#
2200 A6#= .482199394#
2210 A7#=-.193527818#
2220 A8#= .035868343#
2230 X=Z:G#=0
2240 IF X<1 THEN G#=-LOG(X):GOTO 2270
2250 WHILE X>2:X=X-1:G#=G#+LOG(X):WEND
2260 X=X-1
2270 G#=G#+LOG(1+X*(A1#+X*(A2#+X*(A3#+X*(A4#+X*(A5#+X*(A6#+X*(A7#+X*A8#))))))))
2280 L.GAMMA=G#
2290 G=G#
2300 RETURN
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】まず，次元ｎをｎ＝３として，点の総数Ｎを１０→３０→１００としてみる．最近接点の球面距離の平均値は点の密度の増加とともに減少する様子を窺い知ることができる．この結果は直感的にも納得のいくものであろう．
（ｎ＝３）
θ Ｎ＝１０ Ｎ＝３０ Ｎ＝１００
0 0 0 0
5 .0188465 .0554782 .173252
10 .0731504 .203781 .532149
15 .156655 .400173 .817993
20 .260331 .5953 .950972
25 .374052 .754727 .990764
30 .488248 .86596 .998767
35 .595177 .933644 .999882
40 .689595 .970083 .999992
45 .768833 .987641 1
50 .832419 .995291 1
55 .881454 .998332 1
60 .917956 .999447 1
65 .944294 .999827 1
70 .962786 .999948 1
75 .975466 .999985 1
80 .98399 .999996 1
85 .989627 .999999 1
90 .993307 1 1
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【３】次に，点の総数ＮをＮ＝１００として，次元ｎを３→１０→３０としてみる．
（Ｎ＝１００）
θ ｎ＝３ ｎ＝１０ ｎ＝３０
0 0 × ×
5 .173252 × ×
10 .532149 × ×
15 .817993 × ×
20 .950972 × ×
25 .990764 4.46332E-03 ×
30 .998767 .0277094 ×
35 .999882 .0967941 ×
40 .999992 .258247 ×
45 1 .526614 ×
50 1 .807128 8.61132E-03
55 1 .961327 .0369062
60 1 .997164 .188712
65 1 .999943 .590541
70 1 1 .949364
75 1 1 .99966
80 1 1 1
85 1 1 1
90 1 1 1
とくに４よりも大きい次元ｎに関して，θの小さいところで分布関数Ｆの誤差はかなり大きくなった．期待はずれの結果となってしまったが，この原因はプログラムだけにあるのではないと思われる．
コラム「高次元のパラドックスとスターリングの公式（その３）」ではベータ分布がＳn-1球面上で一様に分布する点の配置に密接な関係があることをみてきたが，一様に分布しているといっても大きなｎに対してはほとんどが赤道のかなり近くに位置していることが示された．ここでの計算が示していることも，（直観に反して）大きいｎに対してはＳn-1のほとんどが赤道のかなり近くに位置しているというものであろう．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝