（その２）ではベータ分布がＳn-1球面上で一様に分布する点の配置に密接な関係があることをみてきました．ところが，一様に分布しているといっても大きなｎに対してはほとんどが赤道のかなり近くに位置していることが示されます．今回のコラムでは高次元単位球面上の測度集中について取り上げたいと思います．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】赤道周りへの測度集中
Ｓn-1球面の赤道の周りに幅２ωの帯を設けて，ここに測度の８０％が含まれるようなωを求めてみることにします．すなわち
Ｐ｛ｘ＜Ｓn-1：−ω≦ｘn≦ω｝
＝Ｐ｛ｘ＜Ｓn-1：ｘn^2≦ω^2｝＝０．８
具体的にいうとＳ0＝｛±１｝，Ｓ1は２π（円周）ですから，
ω＝sin(0/8/4･2π)＝０．９５１０５６
ｎ＝３の場合，Ｓ2は４π（表面積）ですが，y=f(x)>0のグラフをｘ軸を中心に回転させてできる曲面の面積は
S[y]=2π∫y(1+(y')^2)^1/2dx
で与えられますから，y=(1-x^2)^(1/2)とおくと
S[y]=2π∫(0,x)dx=2πx=0.8/2･4π
よりω＝０．８と求められます．
ｎ＞３は簡単には求められませんが，（その２）で述べたことにより，x=(x1,x2,･･･,xn)を単位球面Ｓn-1上で一様分布する点とすると，xn^2の分布はベータ分布Beta(1/2,(n-1)/2)となることがわかります．
ベータ分布は不完全ベータ関数と密接に関係していて，その分布関数はガウスの超幾何関数2Ｆ1を使って以下のように表現されます．
Ｆ（ｘ）＝１／ｐｘ^p2Ｆ1(p,1-q,p+1,x)／Β（ｐ，ｑ）
＝１／ｐｘ^p（１−ｘ）^q2Ｆ1(p+q,1,p+1,x)／Β（ｐ，ｑ）
この式からわかるようにベータ分布は区間（０，１）で定義された分布で，（０，１）に制限されているため，ニュートン法などの反復解法を用いてパーセント点を求めるには不向きです．そこで，区間（０，∞）で定義されたＦ分布の上側確率との間には
Ｑ（df1,df2,F)＝Ｉx（df2/2,df1/2,x)，ｘ＝df2/(df2+df1*F)
Ｉx(p,q)＝∫(0,x)t^(p-1)(1-t)^(1-q)dt/B(p,q)
の関係のあることを使って，Ｆ分布に関する計算をしてからベータ分布関数に変換する方法でベータ分布のパーセント点を求めることにします．また，その方が初期値を求めるうえでもを便利です．
ｎ ω
2 .951057
3 .8
4 .687049
5 .6084
6 .550863
7 .506727
8 .471589
9 .442796
10 .418662
20 .291384
30 .236612
40 .204344
50 .182466
60 .166381
70 .153916
80 .143889
90 .135596
100 .12859
この計算が示していることは，（直観に反して）大きいｎに対してはＳn-1のほとんどが赤道のかなり近くに位置しているというものです．後述するように，ｎが大きいときωのオーダーはｎ^(-1/2)になります．
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
1000 DF2=1
1010 FOR I=2 TO 10
1020 P0=.2 :DF1=I-1
1030 GOSUB *F.PERCENT1 :PRINT #3,UUF;
1040 GOSUB *F.PERCENT.HGF:PRINT #3,UUF
1050 X=DF2/(DF2+DF1*UUF) :PRINT #3,I,SQR(X)
1060 NEXT I
1070 DF2=1
1080 FOR I=20 TO 100 STEP 10
1090 P0=.2 :DF1=I-1
1100 GOSUB *F.PERCENT1 :PRINT #3,UUF;
1110 GOSUB *F.PERCENT.HGF:PRINT #3,UUF
1120 X=DF2/(DF2+DF1*UUF) :PRINT #3,I,SQR(X)
1130 NEXT I
1140 END
1150 '
1160 ' ** F-DISTRIBUTION **
1170 '
1180 *F.PERCENT.HGF:
1190 GOSUB *F.PERCENT1:UU=UUF
1200 PROB$="F":GOSUB *NEWTON.HGF
1210 UUF=X
1220 RETURN
1230 '
1240 ' *** ニュートン法 ***
1250 '
1260 *NEWTON.HGF:
1270 EPSL=1E-09
1280 IMAX=50
1290 '
1300 Q0=P0
1310 D0=UU
1320 '
1330 SW=0
1340 ICNT=0
1350 WHILE SW=0
1360 ICNT=ICNT+1:PRINT "*";
1370 IF ICNT>IMAX THEN SW=1
1380 DD=D0:GOSUB 1420:D0=DD
1390 WEND
1400 X=DD
1410 RETURN
1420 '
1430 ' ** 逐次計算 **
1440 '
1450 IF PROB$="F" THEN UUF=D0:GOSUB *F.PROB.HGF
1460 '
1470 F=P0-Q0:G=P1
1480 '
1490 F2=F:D2=D0
1500 IF F1=0 THEN 1530
1510 GOSUB *BISECTION
1520 'GOSUB *DAMPING
1530 '
1540 H=P2:DI=G*G-F*H*2
1550 IF DI<0 THEN DELTA=-F/G:GOTO 1580
1560 IF G<0 THEN DELTA=F*2/(-G+SQR(DI)) ELSE DELTA=F*2/(-G-SQR(DI))
1570 '
1580 DD=D0+DELTA
1590 '
1600 'EPSL=.00001
1610 IF ABS(DD-D0)<EPSL THEN SW=1
1620 'IF ABS((DD-D0)/D0)<EPSL THEN SW=1
1630 F1=F2:D1=D2
1640 RETURN
1650 '
1660 *BISECTION:
1670 IF SGN(F1)=SGN(F2) THEN RETURN
1680 DD=(D1+D2)/2:' bisection
1690 DD=(F1*D2-F2*D1)/(F1-F2):' regula falsi
1700 RETURN 1590
1710 '
1720 *DAMPING:
1730 IF ABS(F1)>ABS(F2) THEN RETURN
1740 DD=D0-DELTA/2:' damping
1750 RETURN 1590
1760 '
1770 ' ** F-DISTRIBUTION **
1780 '
1790 *F.PROB.HGF:
1800 GOSUB *F.PROB.HGF.SUB
1810 P0=1-P0
1820 RETURN
1830 '
1840 ' ** F-DISTRIBUTION **
1850 '
1860 *F.PROB.HGF.SUB:
1870 DFX=DF2/2:DFY=DF1/2
1880 A0=DFX+DFY:B0=1:C0=DFX+1
1890 A=A0:B=B0:C=C0
1900 Z=DF2/(DF2+DF1*UUF):GOSUB *GAUSS.HGF :G1#=G#:' [polynomial expansion]
1910 'Z=DF2/(DF2+DF1*UUF):GOSUB *GAUSS.HGF1:G1#=G#:' [ continued fraction ]
1920 P0=G1#/DFX*Z^DFX*(1-Z)^DFY
1930 Z=DFX+DFY:GOSUB *LOG.GAMMA:G2#=G#
1940 Z=DFX :GOSUB *LOG.GAMMA:G3#=G#
1950 Z=DFY :GOSUB *LOG.GAMMA:G4#=G#
1960 P0=P0*EXP(G2#-G3#-G4#)
1970 'P0=1-P0
1980 '
1990 Z=DF2/(DF2+DF1*UUF)
2000 A=A0+1:B=B0-1:C=C0:G0#=1#
2010 G#=-(G1#*(C0-A0-1)+G0#*(B0-C0))/(A0-B0+1)/(1-Z)
2020 G5#=(G#-G1#)*A0/Z
2030 ZZ=-DF2*DF1/(DF2+DF1*UUF)^2
2040 P10=Z^DFX*(1-Z)^DFY*((DFX/Z-DFY/(1-Z))*G1#+G5#)
2050 P1=-P10*ZZ/DFX*EXP(G2#-G3#-G4#)
2060 '
2070 A=A0+1:B=B0-1:C=C0:G0#=1#
2080 G6#=-G#*A0*(C0-B0*2+(B0-A0-1)*Z)/Z/(1-Z)-G0#*A0*(B0-C0)/Z/(1-Z)-G5#*C0/Z
2090 ZZZ=-DF2*DF1*DF1*2/(DF2+DF1*UUF)^3
2100 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)
2110 P2=-P20/DFX*EXP(G2#-G3#-G4#)
2120 RETURN
2130 '
2140 ' *** 合流型超幾何関数 ***
2150 '
2160 *CONFLUENT.HGF:' [Kummer]
2170 EPSL=1E-09
2180 G0#=0:G#=1
2190 T#=1:K=1
2200 WHILE ABS(G0#-G#)>EPSL
2210 T#=T#*A/C*Z/K
2220 G0#=G#
2230 G#=G#+T#
2240 K=K+1
2250 A=A+1
2260 C=C+1
2270 WEND
2280 G=G#
2290 RETURN
2300 '
2310 ' *** ガウス型超幾何関数 ***
2320 '
2330 *GAUSS.HGF:
2340 EPSL=1E-09
2350 G0#=0:G#=1
2360 T#=1:K=1
2370 WHILE ABS(G0#-G#)>EPSL
2380 T#=T#*A*B/C*Z/K
2390 G0#=G#
2400 G#=G#+T#
2410 K=K+1
2420 A=A+1
2430 B=B+1
2440 C=C+1
2450 WEND
2460 G=G#
2470 RETURN
2480 '
2490 ' *** 合流型超幾何関数の連分数展開 ***
2500 '
2510 *CONFLUENT.HGF1:' 1F1(1,c,x)
2520 NZ=20
2530 G#=0
2540 'FOR ID=NZ TO 1 STEP -1
2550 ID=NZ
2560 WHILE ID>=1
2570 EVEN#=ID/(C+2*ID-1)/(C+2*ID-2)
2580 ODD#=-(C+ID-1)/(C+2*ID)/(C+2*ID-1)
2590 G#=ODD#*Z/(1#+G#)
2600 G#=EVEN#*Z/(1#+G#)
2610 ID=ID-1
2620 WEND
2630 'NEXT ID
2640 G#=-Z/C/(1#+G#)
2650 G#=1#/(1#+G#)
2660 RETURN
2670 '
2680 ' *** ガウス型超幾何関数の連分数展開 ***
2690 '
2700 *GAUSS.HGF1:' 2F1(a,1,c,x)
2710 NZ=20
2720 G#=0
2730 'FOR ID=NZ TO 1 STEP -1
2740 ID=NZ
2750 WHILE ID>=1
2760 ODD#=-(A+ID-1)*(C+ID-2)/(C+2*ID-2)/(C+2*ID-3)
2770 EVEN#=-ID*(C-A+ID-1)/(C+2*ID-1)/(C+2*ID-2)
2780 G#=EVEN#*Z/(1#+G#)
2790 G#=ODD#*Z/(1#+G#)
2800 ID=ID-1
2810 WEND
2820 'NEXT ID
2830 G#=1#/(1#+G#)
2840 RETURN
2850 '
2860 ' *** 対数ガンマ関数 (HASTINGSの多項式近似) ***
2870 '
2880 *LOG.GAMMA:
2890 A1#=-.577191652#
2900 A2#= .988205891#
2910 A3#=-.897056937#
2920 A4#= .918206857#
2930 A5#=-.756704078#
2940 A6#= .482199394#
2950 A7#=-.193527818#
2960 A8#= .035868343#
2970 X=Z:G#=0
2980 IF X<1 THEN G#=-LOG(X):GOTO 3010
2990 WHILE X>2:X=X-1:G#=G#+LOG(X):WEND
3000 X=X-1
3010 G#=G#+LOG(1+X*(A1#+X*(A2#+X*(A3#+X*(A4#+X*(A5#+X*(A6#+X*(A7#+X*A8#))))))))
3020 L.GAMMA=G#
3030 G=G#
3040 RETURN
3050 '
3060 ' ** NORMAL DISTRIBUTION **
3070 '
3080 *NORMAL.PERCENT1:
3090 'P0=1-PP
3100 AA1=2.30753
3110 AA2=.27061
3120 AA3=.99229
3130 AA4=.04481
3140 Q0=.5-ABS(P0-.5)
3150 W=SQR(-2*LOG(Q0))
3160 W1=AA1+AA2*W
3170 W2=1+W*(AA3+W*AA4)
3180 UUN=W-W1/W2
3190 UUN=UUN*SGN(P0-.5)
3200 UU=UUN
3210 RETURN
3220 '
3230 '*NORMAL.PERCENT1:
3240 'P0=1-PP
3250 AA1#=2.515517#
3260 AA2#=.802853#
3270 AA3#=.010328#
3280 AA4#=1.432788#
3290 AA5#=.189269#
3300 AA6#=.001308#
3310 Q0=.5-ABS(P0-.5)
3320 W=SQR(-2*LOG(Q0))
3330 W1=AA1#+W*(AA2#+W*AA3#)
3340 W2=1+W*(AA4#+W*(AA5#+W*AA6#))
3350 UUN=W-W1/W2
3360 UUN=UUN*SGN(P0-.5)
3370 UU=UUN
3380 RETURN
3390 '
3400 ' ** T-DISTRIBUTION **
3410 '
3420 *T.PERCENT1:
3430 GOSUB *NORMAL.PERCENT1
3440 W=UUN*(1+2/(1+8*DF))
3450 UUT=DF*(EXP(W*W/DF)-1)
3460 UUT=SQR(UUT)
3470 UUT=UUT*SGN(P0-.5)
3480 UU=UUT
3490 RETURN
3500 '
3510 ' ** F-DISTRIBUTION **
3520 '
3530 *F.PERCENT1:
3540 IF DF2=1 THEN DF=DF2:P0=1-P0/2:GOSUB *T.PERCENT1 :UUF=1/UUT^2:P0=(1-P0)*2:RETURN
3550 IF DF2=2 THEN UUF=2/DF1/(P0^(-2/DF1)-1):RETURN
3560 FSW=0
3570 IF P0<.5 THEN P0=1-P0:SWAP DF1,DF2:FSW=1
3580 GOSUB *NORMAL.PERCENT1
3590 IF DF2<=3 THEN U=UUN/(DF2^(3/4)) :UUN=UUN*(1.1581-.2296*U+.0042*U*U+.0027*U*U*U)
3600 AA1=2/(9*DF1)
3610 AA2=2/(9*DF2)
3620 W=UUN*UUN
3630 W1=1+AA2*(AA2-W-2)
3640 W2=AA1+AA2-AA1*AA2-1
3650 W3=1+AA1*(AA1-W-2)
3660 W4=SQR(W2*W2-W1*W3)
3670 UUF=(W4-W2)/W1
3680 UUF=UUF^3
3690 IF FSW=1 THEN UUF=1/UUF:P0=1-P0:SWAP DF1,DF2:FSW=0
3700 UU=UUF
3710 RETURN
3720 '
3730 ' ** CHI-SQUARE DISTRIBUTION **
3740 '
3750 *CHI2.PERCENT1:
3760 GOSUB *NORMAL.PERCENT1
3770 AA1=2/(9*DF)
3780 W=1-AA1+UUN*SQR(AA1)
3790 UUX=DF*(W^3)
3800 UU=UUX
3810 RETURN
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】球面上の測度集中
前節では赤道周りへの測度集中について調べてみましたが，ここでは赤道の周りに限定せず，球面Ｓn-1の半分以上を占める任意の部分集合Ａ，すなわち，
Ｐ（Ａ）≧１／２
に対してＡ周りの測度集中を考えます．
ＡtをＡへのユークリッド距離がｔ以下の点ｘ（＜Ｓn-1）の集合とするとき，Ａtでない確率は
１−Ｐ（Ａt）≦２ｅｘｐ（ｰt^2n/2）
であることが示されています．
右辺はｔが増加するにつれて指数関数的に減少するわけですから，したがって，Ａが球面の半分以上を占めるとき，球面上のほとんどすべての点がＡから高々
ｔ＝Ｏ（ｎ^(-1/2)）
の距離しか離れないところにあるというわけです．
［参］Ａが半分以下
Ｐ（Ａ）＝α （０＜α≦１／２）
のとき
１−Ｐ（Ａt）≦２ｅｘｐ（ｰ(t-t0)^2n/2）
が成り立つ．ただし，ｔ0は
２ｅｘｐ（ｰt0^2n/2）＜α
を満たすものとする．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【３】球面上の測度集中の応用
前述したことは，確率論の観点からＳn-1上で定義された確率変数（の和）に関する裾評価を与えてくれます．例としてコルモゴロフ・スミルノフ検定をあげますが，コルモゴロフ・スミルノフ検定はウィルコクソン検定と同じく，分布の形が等しいときに位置に違いがあるかどうかを調べるのに有効なノンパラメトリック法です．この検定には位置が接近していても分布の形が異なれば比較的大きな検出力をもつという特徴があり，分布型の検定にも用いられます．
コルモゴロフ・スミルノフ検定には，見かけは異なるものの実質的には同等ないくつかの表現型がありますが，直観的にわかりやすいのは次のような方法でしょう．それぞれの母集団からの標本を
Ｘ1，Ｘ2，・・・，Ｘn1
Ｙ1，Ｙ2，・・・，Ｙn2
として，このｎ1＋ｎ2個の標本を大きさの順に並べたものを
Ｚ1≦Ｚ2≦・・・≦Ｚ(n1+n2)
とおき，ＺiがＸiのいずれかであるときにはスコアＳi＝１，ＹiのいずれかであるときにはスコアＳi＝−１とします．
スコアの与え方からわかるように，コルモゴロフ・スミルノフ検定はランダム・ウォークなど確率過程と深く関係していて，数学的には面白い性質をもっています．そのため，その解明には難解で高度な数学的知識を必要とします．
たとえば，コルモゴロフ・スミルノフ検定を使うと中央値検定が高い検出力をもつ分布（たとえば，両側指数分布など）に対しては，よい検出力を与えます．しかし，母集団分布が正規分布に従うと仮定した場合，ウィルコクソン検定がｔ検定と同等の検出力をもつための検定効率は０．９５と高いのに対し，正規分布のずれのモデルに関してコルモゴロフ・スミルノフ検定の検定効率は０．６４〜０．７５とかなり小さくなるとされています．（かなり前のことなのですが，コルモゴロフ・スミルノフ検定のｋ標本問題の検定効率に関するDoobの論文を読んでみて，私にはまったく歯がたたなかったことが思い出されました．）
さて，コルモゴロフ・スミルノフ統計量ｋは，確率変数Ｓiの和＝標本の累積度数の差の最大値
ｋ＝max｜ΣＳi｜
と書くことができます．この場合のコルモゴロフ・スミルノフ統計量の漸近分布は計算されていて，
ｖ＝ｋ√（（ｎ1＋ｎ2）／ｎ1ｎ2））
がP{V≧v}となる漸近確率はp=2exp(-2v*v)で与えられます．これが与えられた有意水準より小さくなるか大きくなるかで仮説を棄却・採択します．
この根拠が前節に掲げた
１−Ｐ（Ａt）≦２ｅｘｐ（ｰt^2n/2）
というわけです．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【４】ガウス測度の集中
（その２）では，格子上のランダムウォークについて
(1)ｎステップの後の１次元酔歩が原点からの距離が√ｎのオーダーであること
(2)ｎステップの後のｄ次元酔歩も原点からの距離が√ｎのオーダーであること
が示されましたが，これは拡散現象では一般的にいえることであることがわかります．
格子上のランダムウォークは，ブラウン運動などの拡散モデルとしてよく知られていますが，格子のモデルはブラウン運動の離散化とみなすことができるので，局所的にみると離散的過程であっても，大域的にみると連続空間に分布した連続的なｎ次元ガウス分布とみることができます．
ｎ次元ガウス分布にしたがってランダムに選んだ点の原点からの距離は期待値√ｎの周りに鋭く集中していますから，半径√ｎの球面Ｓn-1上の一様分布からランダムに点を選ぶことに類似しています．このことから，ガウス測度の集中に対する同様な結果については
γ（Ａ）≧１／２
γ（Ａt）≧１−ｅｘｐ（ｰt^2/2）
が成り立つことが知られています．
この不等式には次元ｎが現れないのですが，それは半径√ｎの球面Ｓn-1に対する状況と非常によく似ていることから理解されます．また，係数２は現れませんから，球面の状況と本質的に違うところは係数だけと考えることができます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝