■正規楕円とｔ楕円
コラム「ｎ次元楕円の陰と影」では，Ｘがｎ次元正規分布に従い，その分散共分散行列がΣで表されるとき，ｎ次元楕円
(X-μ)'Σ^(-1)(x-μ)=c^2
の内部の点の全確率ｐが自由度ｎのχ2分布
ｐ＝χ2（ｎ）
で与えられることを示した．
また，それを（Ｘi，Ｘj）平面に投影した場合，全体のｐ％がその内側に入るような楕円を描くには，上記のｃ^2を自由度２のχ2分布（すなわち，指数分布）によって，
c^2=-2ln(p)
と定めればよいこともわかった．
この結果は，確率変数ｘiが標準正規分布Ｎ（０，１）に従うとき，ｘi^2の分布は自由度１のχ2分布，また，ｎ個の変数ｘiがすべてＮ（０，１）に従うならば，Σｘi^2は自由度ｎのχ2分布になる，すなわち，
x〜N(0,1) → x^2〜χ2(1)
xi〜N(0,1) → Σxi^2〜χ2(n)
より，直感的にも理解されるところであろう．
それでは，Ｘがｎ次元ｔ分布に従うとき，当該の問題の解はどのようになるだろうか？
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】ｔ楕円
数学の定理や命題において，２〜３次元で成立することは一般次元でも同じように成立することが多い．しかし，
x〜t(df) → x^2〜{ｔ(df)}^2=F(1,df)
は成り立つものの，これは一般次元では成立しない．決して
xi〜t(df) → Σxi^2〜{ｔ(df)}^(n+1)=F(n+1,df)
などというでたらめを書かないように！．．．というわけで，解を導き出してみたい．
小生が統計の問題を解く際に，いつも頼みの綱にしているいわば命綱が，Kotz,Balakrishnan,Johnson「Continuous univariate distributions 2nd ed.」John Wiley & Sonsである．ところが，どうしたわけか「Continuous multivariate distributions」の新版（2nd ed.）では多変量ｔ分布，Wishart分布がオミットされている．生命線を絶たれたかに思われたのだが，蓑谷千鳳彦「すぐに役立つ統計分布」東京図書を参考にすることができたのは幸運であった．
それによると，自由度ｍのｎ変量ｔ分布の同時確率密度関数は，
f(X,μ,Σ)=Γ((m+n)/2)/Γ(m/2)(mπ)^(n/2)|Σ|^(-1/2)(1+(X-μ)'Σ^(-1)(x-μ)/n)^(-(m+n)/2)
で与えられる．ここで，μは平均ベクトル，Σはｎ×ｎ次の分散共分散行列である．
一般に，ｎ次元楕円
(X-μ)'Σ^(-1)(x-μ)=c^2
を直交座標系（Ｏ：Ｘ1，Ｘ2，Ｘ3，・・・）での関数式で表すと交差項ＸiＸjが出現するため，取り扱いが厄介である．そこで，このｎ次元楕円は，座標変換により，別の直交座標系（ｏ：ｘ1，ｘ2，ｘ3，・・・）において，以下のような標準形
ｘ1^2／ａ^2＋ｘ2^2／ｂ^2＋・・・＋ｘn^2／ｎ^2＝１
で表されるものとする．
正定値対称行列のとき，行列の固有値を
λ1，・・・，λn
とすると，すべての固有値は正であり，
λ1ｘ1^2＋λ2ｘ2^2＋・・・＋λnｘn^2＝ｃ^2 （ｃ^2は定数）
が成り立つから，楕円半径は
ａ^2＝ｃ^2／λ1，ｂ^2＝ｃ^2／λ2，・・・，ｎ^2＝ｃ^2／λn
で表されることになる．
また，このとき，
｜Σ｜＝λ1・・・λn
であり，
｜（Ｘ−μ）’Σ^(-1)（Ｘ−μ）｜
＝｜（Ｘ1−μ1）^2／λ1＋・・・＋（Ｘn−μn）^2／λn｜
＝ｃ^2
の値は不変であるから，変換後の座標で
（Ｘ1−μ1）／√λ1＝ｘ1，・・・，（Ｘn−μn）／√λn＝ｘn
と標準化すると，球面上で一様分布する点となる．
そこで，正規分布の場合と同様に，当該の積分
∫∫・・∫∫ｆ（Ｘ，μ，Σ）ｄＸ1ｄＸ2・・ｄＸn
を直交座標でなく極座標で求めてみることにする．すなわち，被積分関数を原点を中心とする半径ｒの球面上で積分し，次にｒ＝０からｒ＝ｃまで積分すると，半径ｒの球面上で被積分関数は一定値(1+r^2/n)^((m+n)/2)をとり，ｎ次元超球の体積をＶnとすると表面積はｎＶnｒ^(n-1)，また，Ｖn=π^(n/2)/Γ(n/2+1)であるから，
p=∫(0,c)r^(n-1)(1+r^2/n)^((m+n)/2)dr*2/B(n/2,m/2)/m^(n/2)
z=(1+r^2/n)と変数変換すると
p=1/B(n/2,m/2)∫(0,c^2/n)z^(n/2-1)(1+z)^((m+n)/2)dz
このように，球面上で一様分布する点はベータ分布に密接に関係していることが示された．この不完全ベータ関数は自由度（ｎ，ｍ）のＦ分布であり，また，ｃ^2／ｎまでの定積分であるから，全体のｐ％がその内側に入るような楕円を描くには，
ｐ＝ｎＦ（ｎ，ｍ）
より，上記のｃ^2を自由度（ｎ，ｍ）のＦ分布の下側確率×ｎによって定めればよいことが理解されるだろう．
このことは，幾何学的には，
(1)ｔ分布では正規分布よりも一回り大きな正方形を考えなければならない．
(2)正規変量の２乗和を扱うときはχ2円領域，ｔ変量の２乗和を扱うときは，それよりも一回り大きなＦ円領域を考えなければならない．
ということであって，とくに，Ｆ円領域はScheffeの多重比較法の考え方に類似しているのである．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】ｎ次元正規楕円と相関を考慮したモンテカルロ法
正規乱数を発生させる方法に，ボックス・ミューラー（Box-Muller）法がある．ボックス・ミューラー法は，レイリー分布を応用して正規乱数を発生させる方法であるから，同じ標的問題の解であるマクスウェル分布など一般化したカイ分布（カイ２乗分布の平方根の分布）を用いることによって，ｎ個の一様乱数からｎ個の互いに独立な正規乱数を容易に作り出すことができる．
また，そのように凝った手を用いずとも，ｎ個の正規乱数を単純に組み合わせるだけであっても，そのランダム性は失われないわけであるから，相関のないｎ次元正規乱数を発生させることは可能である．
しかし，通常，変数同士には相関があり，どうしても相関を考慮に入れたモンテカルロ法が必要になる場合がある．そこで，この節では分散共分散行列のコレスキー分解に基づいた多変量正規乱数発生法について紹介することにしたい．
ｎ次元正規分布の同時確率密度関数は
f(X,μ,Σ)=(2π)^(n/2)Σ^(-1/2)exp(-(X-μ)'Σ^(-1)(x-μ)/2}
と表すことができるが，Σが分散共分散行列であり，(X-μ)'Σ^(-1)(x-μ)=c^2が等確率楕円である．
この等確率楕円が描けるような多変量正規乱数発生法の原理は簡単で，互いに独立な標準正規乱数Ｘ1，Ｘ2，・・・，Ｘnから，
Ｙ1＝μ1＋ａ11Ｘ1
Ｙ2＝μ2＋ａ21Ｘ1＋ａ22Ｘ2
・・・・・・・・・・・・・・
Ｙn＝μn＋ａn1Ｘ1＋ａn2Ｘ2＋・・・＋ａnnＸn
とおくと，Ｙiは正規変数の線形結合であるから平均値μiの正規分布することは明らかだし，ＹiとＹjが相関することも理解できるだろう．
その際の分散・共分散σijは，
σij＝ａi1ａj1＋ａi2ａj2＋・・・＋ａijａjj （ｉ≧ｊ）
に等しい．したがって，あとはａijが求められればよいことになるが，それには正値対称行列Σのコレスキー分解
Σ＝ＬＬ’ （Ｌは下三角行列，Ｌ’はＬの転置行列）
によって
Ｙ＝ＬＸ＋μ
と置けばよいことが示される．詳細については，宮武・脇本共著「乱数とモンテカルロ法」共立出版を参照されたい．
下三角行列とは行列要素の上方の三角部が０で対角要素と下方の三角部だけからなる行列であるが，最小２乗法のような正値対称行列ではＬＵ分解やコレスキー分解で効率的に連立１次方程式を解くことができ，計算スピードも加速される．また，修正コレスキー分解とは
Σ＝ＬＤＬ’ （Ｄは対角行列）
のように三角化と対角化を行うマトリックス計算法である．
ガウス消去法やガウス・ジョルダン法は，丸め誤差がいたずらをするので基本的には不安定であるが，コレスキー分解は計算精度の点でも優れていることが知られている．コレスキー分解では手間もガウス消去法の半分で済むから，格段と能率がよい解法なのであるが，相関を考慮に入れた多変量正規乱数発生にも応用できるのである．意外な使い道であり，予期せざる効能といってもよいだろう．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
