ｎ次元楕円の投影問題は，元々，測定データのなかに飛び離れた値（外れ値：outlier）があるとき，棄却すべきか採択すべきかを判断する統計検定に由来しています．それを視覚化したものが確率楕円なのですが，確率楕円はたとえば９５％のデータ点を含むような楕円として描かれます．このような確率構造が入るのは，棄却検定という統計学上の問題が設定されているからなのです．

　

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

　　(1)確率構造の入った超楕円体の投影問題を扱ってみること

にしました．また，超楕円体と超直方体の比較を行うために，

　　(2)確率構造の入った超直方体の投影問題

も扱っています．

　

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

　

【１】集中楕円（concentrative ellipsoid）

　

　主成分分析（ＰＣＡ）における第１主成分を表す直線（直交回帰直線）は，データ点のそれぞれからその直線に下した垂線も長さの２乗和を最小となるようにしたものです．これは第１主成分の分散を最大にすることに対応しています．直交回帰直線は算術平均を通り，計測値の集団の確率楕円の長軸に一致しますから，算術平均を中心として直交回帰直線を長軸とした楕円を描くと確率楕円が得られます．

　

　主成分分析は多変量データ解析において最も基本的な方法の一つですが，理解のためには固有値や固有ベクトルの知識が必要になります．これらについては（その１）をご参照願います．

　

［１］ｎ次元正規分布の集中楕円

　

　ｎ次元正規分布において，ｎ個の確率変数（ｘ1，・・・，ｘn）の同時確率は

　　ｆ（ｘ，μ，Σ）＝（２π）^(n/2)Σ^(-1/2)ｅｘｐ｛−（ｘ−μ）’Σ^(-1)（ｘ−μ）／２｝

によって与えられます．ここで，ｘ＝（ｘ1，・・・，ｘn）’はデータベクトル，μ＝（μ1，・・・，μn）’は平均ベクトル，Σはｎ×ｎ次の分散共分散行列を表します．

　

　この節では，ｎ次元空間内の点

　　（μ1，・・・，μn）

を中心として，ｎ次元正規分布するデータ点に対して，全体のｐ％がその内側に入るような確率のｎ次元楕円を定めるという問題を考えることにします．

　

　上式の指数部の引数を−ｃ^2／２とおいてみると，楕円面上の点はすべて

　　ｆ（ｘ，μ，Σ）＝（２π）^(n/2)Σ^(-1/2)ｅｘｐ｛−ｃ^2／２｝

という確率をもつ等確率楕円

　　（ｘ−μ）’Σ^(-1)（ｘ−μ）＝ｃ^2

を表すことになります．

　

　また，（その１）のｘ’Ｈｘ＝１との対応でいえば，

　　ｘ　←→　（ｘ−μ）

　　Ｈ　←→　Σ^(-1)

　　１　←→　ｃ^2

ということになりますが，統計学では，ｘ’Ｈｘよりも分散・共分散行列Σを用いたｘ’Σ^(-1)ｘの形式の方が愛されているようです．以下，ＨとΣを適宜使い分けることにします．

　

　この等確率楕円面の内部の点の全確率は，多重積分

　　∫∫・・∫∫ｆ（ｘ，μ，Σ）ｄｘ1ｄｘ2・・ｄｘn

により計算できます．とはいってもこのままでは計算は困難です．一般に，ｎ次元楕円

　　(x-μ)'Σ^(-1)(x-μ)=c^2

を直交座標系（Ｏ：ｘ1，ｘ2，ｘ3，・・・）での関数式で表すと交差項ｘiｘjが出現するため，取り扱いが厄介だからです．

　

　そこで，このｎ次元楕円は，座標変換により，別の直交座標系（ｏ：Ｘ1，Ｘ2，Ｘ3，・・・）において，以下のような標準形

　　Ｘ1^2／ａ^2＋Ｘ2^2／ｂ^2＋・・・＋Ｘn^2／ｎ^2＝１

で表されるものとします．ここで，Ｘi＝ｘi−μi，また，ａ，ｂ，・・・は楕円半径を表します．

　

　分散・共分散行列Σは正定値対称行列ですから，その固有値を

　　λ1，・・・，λn

とすると，すべての固有値は正であり，

　　Ｘ1^2／λ1＋Ｘ2^2／λ2＋・・・＋Ｘn^2／λn＝ｃ^2　　　（ｃ^2は定数）

が成り立ちますから，楕円半径は

　　ａ^2＝λ1ｃ^2，ｂ^2＝λ2ｃ^2，・・・，ｎ^2＝λnｃ^2

で表されることになります．

　

　また，このとき，

　　｜Σ｜＝λ1・・・λn

であり，

　　｜（ｘ−μ）’Σ^(-1)（ｘ−μ）｜

　＝｜Ｘ1^2／λ1＋・・・＋Ｘn^2／λn｜

　＝ｃ^2

の値は不変ですから，変換後の座標で

　　（ｘ1−μ1）／√２λ1＝ｙ1，・・・，（ｘn−μn）／√２λn＝ｙn

と標準化すると，（ｙ1，・・・，ｙn）は球面上で一様分布する点，一方，当該の多重積分は

　　∫∫・・∫∫ｆ（ｘ，μ，Σ）ｄｘ1ｄｘ2・・ｄｘn

　＝∫∫・・∫∫π^(-n/2)ｅｘｐ（−ｙ1^2−・・・−Ｘn^2）ｄｙ1ｄｙ2・・ｄｙn

となります．

　

　　Ｉ＝∫・・∫ｅｘｐ（−ｙ1^2−・・・−ｙn^2）ｄｙ1ｄｙ2・・ｄｙn

とおくと，この積分は多重ガウス積分であり，

　　Ｉ＝｛∫ｅｘｐ（−ｙ^2）ｄｙ｝^n

　

　ここで，直交座標でなく極座標，すなわち，被積分関数を原点を中心とする半径ｒの球面上で積分し，次にｒ＝０からｒ＝ｃまで積分をするのです．すると，半径ｒの球面上で被積分関数は一定値exp(-r^2)をとり，ｎ次元超球の体積をＶnとすると表面積はｎＶnｒ^(n-1)ですから，

　　Ｉ＝∫(0,c)ｅｘｐ（−ｒ^2）ｎＶnｒ^(n-1)ｄｒ

　　　＝ｎＶn∫(0,c)ｒ^(n-1)ｅｘｐ（−ｒ^2）ｄｒ

　

　ｚ＝ｒ^2と置換すると，ｄｚ＝２ｒｄｒより，

　　Ｉ＝ｎＶn／２∫(0,c^2)ｚ^(n/2-1)ｅｘｐ（−ｚ）ｄｚ

また，

　　Ｖn=π^(n/2)/Γ(n/2+1)

と表されますから，

　　Ｉ＝π^(n/2)γ(n/2,c^2)／Γ(n/2)

　

　このように不完全ガンマ関数／ガンマ関数が現れましたから，結局，求めるｃ^2はガンマ分布の分布関数が

　　γ(n/2,c^2)／Γ(n/2)

と表されることより，自由度ｎのχ^2分布の下側ｐ％点で与えられることになります．

　　ｃ^2＝χ^2（ｎ：ｐ）

　

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

　

［２］ｎ次元ｔ分布の集中楕円

　

　上の結果は，確率変数ｘが標準正規分布Ｎ（０，１）に従うとき，ｘ^2の分布は自由度１のχ^2分布，また，ｎ個の変数ｘiがすべてＮ（０，１）に従うならば，Σｘi^2は自由度ｎのχ^2分布になる，すなわち，

　　x〜N(0,1) → x^2〜χ^2(1)

　　xi〜N(0,1) → Σxi^2〜χ^2(n)

より，直感的にも理解されるところです．

　

　それでは，ｘがｎ次元ｔ分布に従うとき，当該の問題の解はどのようになるのでしょうか？　確率変数ｘが自由度ｄｆのｔ分布に従うとき，ｘ^2の分布は自由度（１，ｄｆ）のＦ分布となります．

　　x〜t(df) → x^2〜{ｔ(df)}^2=F(1,df)

　

　数学の定理や命題において，１〜３次元で成立することは一般次元でも同じように成立することが多いのですが，

　　x〜t(df) → x^2〜{ｔ(df)}^2=F(1,df)

は成り立つものの，これは一般次元では成立しません．決して

　　xi〜t(df) → Σxi^2〜{ｔ(df)}^(n+1)=F(n,df)

などというでたらめを書かないように！．．．というわけで，解を導き出してみることにします．

　

　自由度ｄｆのｎ変量ｔ分布の同時確率密度関数は，

　　f(x,μ,Σ)=Γ((df+n)/2)/Γ(df/2)(dfπ)^(n/2)|Σ|^(-1/2)(1+(x-μ)'Σ^(-1)(x-μ)/n)^(-(df+n)/2)

で与えられます．

　

　変換後の座標で

　　（ｘ1−μ1）／√λ1＝ｙ1，・・・，（ｘn−μn）／√λn＝ｙn

と標準化して，正規分布の場合と同様に求めてみると，半径ｒの球面上で被積分関数は一定値(1+r^2/n)^((df+n)/2)をとりますから，

　　p=∫(0,c)r^(n-1)(1+r^2/n)^((df+n)/2)dr*2/B(n/2,df/2)/df^(n/2)

z=(1+r^2/n)と変数変換すると

　　p=1/B(n/2,df/2)∫(0,c^2/n)z^(n/2-1)(1+z)^((df+n)/2)dz

　

　このように，球面上で一様分布する点はベータ分布に密接に関係していることが示されます．この不完全ベータ関数は自由度（ｎ，ｄｆ）のＦ分布であり，また，ｃ^2／ｎまでの定積分ですから，全体のｐ％がその内側に入るような楕円を描くには，ｐを下側確率として

　　ｃ^2／ｎ＝Ｆ（ｎ，ｄｆ：ｐ）

したがって，上記のｃ^2を自由度（ｎ，ｄｆ）のＦ分布の下側ｐ％点×ｎによって定めればよいことが理解されます．

　　xi〜t(df) → Σxi^2〜nF(n,df)

　

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

　

［３］ｎ次元集中楕円のｍ次元空間への投影

　

　ｘがｎ次元正規分布に従い，その分散共分散行列がΣで表されるとき，ｎ次元楕円

　　(x-μ)'Σ^(-1)(x-μ)=c^2

の内部の点の全確率ｐが自由度ｎのχ2分布

　　ｃ^2＝χ^2（ｎ：ｐ）

また，ｘがｎ次元ｔ分布に従うときは，

　　ｃ^2＝ｎＦ（ｎ，ｄｆ：ｐ）

で与えられることを示しました．

　

　以上のことは，

　　(1)ｔ分布では正規分布よりも一回り大きな領域を考えなければならない．

　　(2)正規変量の２乗和を扱うときはχ^2円領域，ｔ変量の２乗和を扱うときは，それよりも一回り大きなＦ円領域を考えなければならない．

ということを意味しています．

　

　次に，これらの結果をｍ次元空間に投影することを考えましょう．

　　ｙ1＝［ｘ1−μ1，・・・，ｘm−μm］

　　ｙ2＝［ｘm+1−μm+1，・・・，ｘn−μn］

の同時分布：ｆ（ｙ1，ｙ2）を考えることになるわけですが，

　　Ｈ＝Σ^(-1) ＝［Ａ，Ｂ］　　　Σ＝［Ｏ，Ｐ］

　　　　　　　　　［Ｃ，Ｄ］　　　　　［Ｑ，Ｒ］

また，

　　Ｖ＝Ｄ−ＣＡ^(-1)Ｂ

と書くことにすると，

　　｜Ｈ｜＝｜Ａ｜｜Ｄ−ＣＡ^(-1)Ｂ｜＝｜Ａ｜｜Ｖ｜

が得られます．

　

　これより，

　　（ｘ−μ）’Σ^(-1)（ｘ−μ）＝ｙ1’Ａｙ1＋ｙ2’Ｖｙ2

　　ｆ（ｙ1，ｙ2）＝｜Ａ｜^(1/2)／（２π）^(m/2)ｅｘｐ（−ｙ1’Ａｙ1／２）×｜Ｖ｜^(1/2)／（２π）^((n-m)/2)ｅｘｐ（−ｙ2’Ｖｙ2／２）

すなわち，ｙ1，ｙ2は互いに独立な多変量正規分布に従うことが理解されます．

　

　このことは直観的にも明らかですから，答えはもうおわかりでしょうが，全体のｐ％がその内側に入るようなｎ次元の集中楕円をｍ次元空間に投影する場合，正規分布では

　　ｃ^2＝χ^2（ｍ：ｐ）

ｔ分布では，

　　ｃ^2＝ｍＦ（ｍ，ｄｆ：ｐ）

と定めればよいことがわかります．

　

　また，このとき，ｍ≦ｎですから，

　　χ^2（ｍ：ｐ）≦χ^2（ｎ：ｐ）

　　ｍＦ（ｍ，ｄｆ：ｐ）≦ｎＦ（ｎ，ｄｆ：ｐ）

が成り立ちます．

　

　例をあげて説明しましょう．ｎ次元正規分布に従う確率変数ｘをｎ次元楕円を（ｘi，ｘj）平面に投影する場合，当該の楕円を描くには，ｃ^2を自由度２のχ^2分布（すなわち，指数分布）によって，

　　ｃ^2＝−２ｌｎ（１−ｐ）

と定めます．例えば，下側確率ｐ＝０．９５，０．９９を代入すると，それぞれｃ^2＝６．０，９．２．

　

　逆に，ｃ＝２のとき，ｐ＝０．８６５となりますから，λ1，λ2を小行列Ｏの固有値として，軸の長さが２√λ1，２√λ2の楕円を描くと，内側に８６．５％が含まれることになるのです．

　

　なお，これまでの説明で，（ｎ次元）正規分布やｔ分布を仮定することが難しいのではと思われた読者もおられるかもしれません．しかし，集中楕円では，厳密な分布にこだわる必要はなく，平均値付近にデータが集中しほぼ左右対称になるような場合には正規分布やｔ分布で近似してもそれほど違いを生じません．

　

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

　

　以下に２次元・３次元の場合の当該の値を掲げますが，ｍ次元空間に投影するとはいっても，最終的に表示するのはモニタ画面上ですから平面です．したがって，必要とされるのはｍ＝２のときの値ということになります．

　

［４］２次元平面への投影（ｍ＝２）

　

ａ）ｃ^2値

　　　　　　　　　　p=0.95　　　　　　 p=0.99

正規分布　　　　　　5.99146　　　　　　 9.21034

ｔ分布(df=10)　　　 8.20564　　　　　　15.1189

ｔ分布(df=20)　　　 6.98566　　　　　　11.6979

ｔ分布(df=30)　　　 6.63166　　　　　　10.7807

　

ｂ）ｐ値

　　　　　　　　　　ｃ＝１　　　　　　ｃ＝２　　　　　　ｃ＝３

正規分布　　　　　　.393469　　　　　 .864665　　　　　 .988891

ｔ分布(df=10)　　　 .379079　　　　　 .814066　　　　　 .959614

ｔ分布(df=20)　　　 .386086　　　　　 .838494　　　　　 .97566

ｔ分布(df=30)　　　 .388504　　　　　 .84702　　　　　　.980463

　

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

　

［５］３次元空間への投影（ｍ＝３）

　

ａ）ｃ^2値

　　　　　　　　　　p=0.95　　　　　　 p=0.99

正規分布　　　　　　7.77862　　　　　 11.3691

ｔ分布(df=10)　　　11.1161　　　　　　19.7827

ｔ分布(df=20)　　　 9.26607　　　　　 14.827

ｔ分布(df=30)　　　 8.73465　　　　　 13.5393

　

ｂ）ｐ値

　　　　　　　　　　ｃ＝１　　　　　　ｃ＝２　　　　　　ｃ＝３

正規分布　　　　　　.19653　　　　　　.739686　　　　　 .971316

ｔ分布(df=10)　　　 .19692　　　　　　.682058　　　　　 .918658

ｔ分布(df=20)　　　 .196686　　　　　 .708782　　　　　 .945894

ｔ分布(df=30)　　　 .196546　　　　　 .718561　　　　　 .954725

　

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

　

【２】集中楕円の計算プログラム

　

　この節では，前節［４］，［５］に掲げた

　　(1)下側確率ｐ→パーセント点ｃ^2

　　(2)パーセント点ｃ^2→下側確率ｐ

の近似計算を行うプログラムを紹介します．

　

　近似計算とはいっても，ｎ次元正規分布に従う場合の平面投影（ｍ＝２）のときは正確な値が求められますし，ｍ≧３の場合であっても相対誤差は％オーダーですから，実用的には十分な有効数字があります．

　

1000 M=2

1010 P0=.95:PP=1-P0

1020 DF=M:GOSUB *CHI2.PERCENT:PRINT UUX

1030 '

1040 DF=M

1050 FOR C=1 TO 3

1060　 UUX=C*C:GOSUB *CHI2.PROB:PRINT 1-PP

1070 NEXT C

1080 '

1090 M=2

1100 P0=.95:PP=1-P0:DF2=10

1110 DF1=M:GOSUB *F.PERCENT:PRINT UUF*DF1

1120 '

1130 DF1=M

1140 FOR C=1 TO 3

1150　 UUF=C*C/DF1:GOSUB *F.PROB:PRINT 1-PP

1160 NEXT C

1170 END

1180 '

1190 ' ** NORMAL DISTRIBUTION **

1200 '

1210 *NORMAL.PERCENT:

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

1230 UUN=SQR(XXX*(2.06118-5.72622/(XXX+11.6406)))

1240 IF PP>.5 THEN UUN=-UUN

1250 RETURN

1260 '

1270 ' ** T-DISTRIBUTION **

1280 '

1290 *T.PERCENT:

1300 GOSUB *NORMAL.PERCENT

1310 UUN=ABS(UUN)

1320 UU2=UUN*UUN

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

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

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

1360 UUT=UUN*(1+AA1+AA2+AA3)

1370 IF PP>.5 THEN UUT=-UUT

1380 RETURN

1390 '

1400 ' ** F-DISTRIBUTION **

1410 '

1420 *F.PERCENT:

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

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

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

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

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

1480 GOSUB *NORMAL.PERCENT

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

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

1510 RETURN

1520 '

1530 ' ** CHI-SQUARE DISTRIBUTION **

1540 '

1550 *CHI2.PERCENT:

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

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

1580 GOSUB *NORMAL.PERCENT

1590 UUX=2/(9*DF)

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

1610 UUX=UUX^3*DF

1620 RETURN

1630 '

1640 ' ** NORMAL DISTRIBUTION **

1650 '

1660 *NORMAL.PROB:

1670 AA1=.196854

1680 AA2=.115194

1690 AA3=.000344

1700 AA4=.019527

1710 W=ABS(UUN)

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

1730 P0=P0^4

1740 P0=1-.5/P0

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

1760 PP=1-P0

1770 RETURN

1780 '

1790 ' ** T-DISTRIBUTION **

1800 '

1810 *T.PROB:

1820 W1=(1-2/(3+8*DF))

1830 W2=DF*LOG(UUT*UUT/DF+1)

1840 UUN=W1*SQR(W2)

1850 GOSUB *NORMAL.PROB

1860 PP=1-P0

1870 RETURN

1880 '

1890 ' ** F-DISTRIBUTION **

1900 '

1910 *F.PROB:

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

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

1940 FSW=0

1950 IF UUF=0 THEN P0=0:GOTO 2060

1960 IF UUF<1 THEN UUF=1/UUF:SWAP DF1,DF2:FSW=1

1970 AA1=2/(9*DF1)

1980 AA2=2/(9*DF2)

1990 W=UUF^(1/3)

2000 W1=W+AA1-W*AA2-1

2010 W2=AA2*W*W+AA1

2020 UUN=W1/SQR(W2)

2030 IF DF2<=3 THEN UUN=UUN*(1+.08*(UUN^4)/(DF2^3))

2040 GOSUB *NORMAL.PROB

2050 IF FSW=1 THEN P0=1-P0:SWAP DF1,DF2:FSW=0

2060 PP=1-P0

2070 RETURN

2080 '

2090 ' ** CHI-SQUARE DISTRIBUTION **

2100 '

2110 *CHI2.PROB:

2120 IF DF=2 THEN PP=EXP(-UUX/2):RETURN

2130 AA1=2/(9*DF)

2140 AA2=UUX/DF

2150 W=AA2^(1/3)

2160 W1=W+AA1-1

2170 UUN=W1/SQR(AA1)

2180 GOSUB *NORMAL.PROB

2190 PP=1-P0

2200 RETURN

　

　近似計算ではどうしても満足できず，正確な値を保証したい方は，コラム「超幾何関数を用いた確率分布の計算」に掲げたプログラムをご利用下さい．以下にその主要部分だけ掲げておきます．

　

1010 FOR I=1 TO 20

1020 DF=I

1030 PRINT:PRINT DF

1040 C=1:UUX=C*C:GOSUB *CHI2.PROB.HGF:PRINT P0

1050 C=2:UUX=C*C:GOSUB *CHI2.PROB.HGF:PRINT P0

1060 C=3:UUX=C*C:GOSUB *CHI2.PROB.HGF:PRINT P0

1110 NEXT I

1120 END

1130 '

1140 ' ** CHI-SQUARE DISTRIBUTION **

1150 '

1160 *CHI2.PROB.HGF:

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

1180 A=A0:C=C0

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

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

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

1220 Z=UUX/2

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

1240 P0=LOG(G1#)-G2#+DF/2*LOG(Z)-UUX/2

1250 P0=EXP(P0)

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

　

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

　

【３】集中多面体（concentrative polyhedron）

　

　まず，復習から始めたいと思います．【１】【２】では集中楕円の問題を扱いましたが，多重ガウス積分は不完全ガンマ関数とになることより，ｎ次元超球では原点からの距離がｃの超球内に納まる確率は自由度ｎのχ^2分布で与えられることが計算されました．

　

　また，２次元平面に投影した場合，自由度２のχ^2分布はすなわち指数分布ですから，全体のｐ％がその内側に入るような超球の半径ｃは

　　ｃ^2＝−２ｌｎ（１−ｐ）

によって与えられることがわかりました．例えば，ｃ＝１のときｐ＝０．３９４，ｃ＝２のときｐ＝０．８６５となり，半径が２の円を描くとその内側には８６．５％が含まれることになります．

　

　これらは，それぞれ１辺が２，４の正方形の内側に入る確率，

　　ｐ＝（０．６８３）^2＝０．４６６＞０．３９４

　　ｐ＝（０．９５４）^2＝０．９１０＞０．８６５

よりもさらに小さくなります．

　

　ここで，０．６８３，０．９５４はそれぞれ（１次元）正規分布の１σ，２σ点における中心確率です．すなわち，正規分布では，区間

　　［μ−σ，μ＋σ］に６８．３％，

　　［μ−２σ，μ＋２σ］に９５．４％，

　　［μ−３σ，μ＋３σ］に９９．７％

の観測値が入ります．３σ点における中心確率は０．９９７となりますから，ほとんどの観測値が［μ−３σ，μ＋３σ］に入ることを利用して，工場では品質管理を行っています．それが３σ法で，有用なＱＣテクノロジーの１つになっています．

　

　一方，ｄ次元正規分布は，

　　p(x1,x2,x3,･･･,xd)=1/(2πσ^2)^(d/2)exp{-(x1^2+x^2+･･･+xd^2)/2σ^2}

で与えられますが，多次元正規分布の場合，低次元の場合とは対照的に密度の裾にあたる領域に大部分のデータが存在することになります．

　

　１次元であれば３σを外れる確率はわずか千に三つにすぎないのですが，高次元化に伴い次第に中心部は過疎化し，１０次元では半分以上が３σの郊外に移り住み，３０次元では９９％以上が裾の領域に集まるという結果になるのです．

　

　さて，変数ｘiはそれぞれ独立に（１次元）正規分布またはｔ分布に従うと考えられますから，この節では全体のｐ％がその内側に入るような超立方体を描く方法について考えてみることにします．

　

　この問題は，統計学ではよく知られたボンフェローニ・シダック問題であって，１回の検定の有意水準をｐ0とすると，ｎ回の検定を繰り返すことによって生じる推論全体での有意水準をｐに保つには，

　　ｐ＝１−（１−ｐ0）^n

より

　　ｐ0 ＝１−（１−ｐ）^(1/n)　　　（シダック解）　

　　　　≒ｐ／ｎ　　　　　　　　　　（ボンフェローニ解）

と求められることと，全く同じ解になっています．

　

　すなわち，ｐ0 ＝１−（１−ｐ）^(1/n)として，ｃを正規分布のｐ0パーセント点

　　ｃ＝Ｎ（ｐ0）

あるいは，ｔ分布のｐ0パーセント点

　　ｃ＝ｔ（ｄｆ：ｐ0）

に定めてから，２次元平面上に投影すればよいことになるのですが，【４】のサブルーチンを使うならば，

　　PP=1-(1-PP)^(1/M):GOSUB *NORMAL.PERCENT:PRINT UUN

といった書式になります．集中楕円の描き方とは，かなり異なることがおわかり頂けたでしょうか？　

　

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

　

　誤解を招くといけませんから，ニワトリが先かタマゴが先かという話をしておきます．統計学では，個々の母数に対する個別の信頼区間は容易に作ることができますが，たとえば，２つ（３つ）の母数が存在する場合，その同時信頼区間は長方形（直方体）領域で与えられるものではなく，楕円（楕円体）領域となります．

　

　すなわち，統計学では，集中楕円体として得られるはずの母数空間を集中多面体で近似することが頻繁に行われているのです．厳密な意味では，集中楕円体が基にあり，集中多面体は集中楕円体を近似するための単なる思いつきの産物ということになりましょう．

　

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