ｔ分布の乱数発生アルゴリズム

　点推定よりも区間推定できるほうが「精度保証つき」計算という面からいって格段と優れていることはいまさらいうまでもありません．

　最小２乗法において，パラメータの信頼区間は分散・共分散行列の対角要素から求めることができます．パラメータの誤差の導出はやや難解ですので，おしつけがましく結果だけを示しておきますが，パラメータθjの１００（１－α）％信頼区間は，

　ここで，Δは分散・共分散行列と呼ばれるｍ次の対称行列で，ｇを行ベクトルとするとき，ｇΔｇ’が正定値２次形式になるという性質をもっています．また，最小２乗推定量θjは最良線形不偏推定量（ＢＬＵＥ），すなわち，線形不偏推定量の中で分散が最小のもの（最小分散不偏推定量）になっています．これらに関する参考文献は多数ありますから，天下り式で我慢できないかたはぜひ参照してください．

　ところで，最小２乗法によって推定されるパラメータの値を使ってモンテカルロ・シミュレーションを行う必要に迫られるときがあります．その際，パラメータの分布はｔ分布に従うわけですから，正規乱数ではなくて，ｔ乱数を発生させなければなりません．

　しかし，ｔ乱数がモンテカルロ・シミュレーションに用いられることはほとんどなく，通常は正規乱数で代用されています．モンテカルロ法の精度はラフなので厳密さにこだわる必要はなかろうというのがおおかたの意見かもしれませんが，今回のコラムでは，モンテカルロ法が粗雑であるからこそ厳密にという考えから，ｔ乱数を用いたモンテカルロ法の話題を取り上げたいと思います．

　ｔ乱数を発生させるには，正規乱数とカイ２乗乱数を組み合わせてつくるのがもっとも近道と思われます．すなわち，ｘが標準正規分布に従い，ｙが自由度dfのχ^2分布に従うとき，

　uを区間[0,1]の一様乱数とすると，それから正規乱数zを発生させる方法はいくつかあって，たとえば，

とすると平均値θ，自由度ｋのガンマ乱数が得られます．なお，ｋが大きいときにはこの方法は効率的ではなく，またｋが整数でないときには適用できません．

　したがって，自由度dfの奇偶によって，df/2が半整数（整数+1/2）あるいは整数のガンマ乱数を作ればよいことになります．自由度が１のガンマ分布は平均１の指数分布ですが，自由度1/2は標準正規分布をする確率変数をｚとするとき，ｚ^2／２の分布に対応します．

　ここで，ガンマ分布の再生性「変数ｘiがガンマ分布G(αi，β)にしたがうとき，Σxiはガンマ分布G(Σαi，β)に従う．」すなわち，ガンマ分布にしたがう変数の和もガンマ分布になるという性質をもっていますから，カイ２乗乱数は指数乱数および正規乱数を使って容易に作れることになります．

　(1)～(3)により，自由度dfのｔ分布が得られたことになりますが，この分布の平均は０，分散はdf/(df-2)で与えられます．

と変数変換すると，平均０，標準偏差１の分布が得られます．これを規格化（standardizationあるいはnormalization）といいます．ｔ分布の場合，

　逆に，もしｔが標準ｔ乱数すなわち平均０，標準偏差１のｔ乱数ならば，平均ｍ，標準偏差ｓのｔ乱数を発生させるには，

　2310～2500行では，正規乱数とガンマ乱数を組み合わせてｔ乱数を発生させています．(1)～(4)の考え方は20行ほどのプログラムとしてまとめることができました．

　たとえば，パラメータ数が３の場合，通常行われている（相関を考慮しない）モンテカルロ法では，３個の標準正規乱数z1,z2,z3を用いて，

　しかし，通常，パラメータ同士には相関があり，どうしても相関を考慮に入れたモンテカルロ法が必要になる場合があります．分散・共分散行列はすべてのパラメータが独立と考えられる場合は対角行列，また，パラメータ間に強い相関がある場合には非対角要素は無視できなくなり対称行列となるからです．

　コラム「正規楕円とｔ楕円」では，分散・共分散行列のコレスキー分解に基づいた多変量正規乱数発生法について紹介しましたが，しばらくの間，それを再録することにします．

　この等確率楕円が描けるような多変量正規乱数発生法の原理は簡単で，互いに独立な標準正規乱数Ｘ1，Ｘ2，・・・，Ｘnから，

とおくと，Ｙiは正規変数の線形結合であるから平均値μiの正規分布することは明らかだし，ＹiとＹjが相関することも理解できるだろう．

に等しい．したがって，あとはａijが求められればよいことになるが，それには正値対称行列Σのコレスキー分解

と置けばよいことが示される．詳細については，宮武・脇本共著「乱数とモンテカルロ法」共立出版を参照されたい．』

　すなわち，相関を考慮に入れたモンテカルロ法でも，必要となる乱数の数は同じ３個であって，唯一異なる点は，

が正規分布することは明らかなのに対し，ｔ分布の場合，再生性は成り立たないわけですから，ｔ分布の線形結合はｔ分布とはなりません．これ以外の議論については正規分布であってもｔ分布であっても同じですから，唯一，この点だけが引っかかってしまいました．

　すなわち，複数の独立なｔ分布に従う乱数の線形結合を使って，多変量のものをつくるというわけにはいかないのではと思われます．とはいっても，具体的にどうすれば多変量ｔ乱数を発生できるかについては思い当たりませんでした．

　しかし，ａn1Ｘ1＋ａn2Ｘ2＋・・・＋ａnnＸnは名目的には線形結合であっても，実質的には乱数１個分に相当するわけですから，近似的な多変量ｔ分布が得られるはずです．この点については，厳密に証明したわけではありませんが，そのように考えないと整合性が保たれないのです．したがって，相関を考慮した多変量ｔ分布の乱数発生でも，ｔを標準ｔ乱数とすると，

　私はこれまで多変量のｔ分布というものを試したことがないので詳しいことは何一つわからないし，確たる自信もないのですが，擬似的な多変量ｔ乱数ということで，ひとつご勘弁を願います．

　ｔ乱数を発生させる基本的な部分はすでにできているので，あとは分散・共分散行列を入力して，コレスキー分解する部分をつくれば完成です．そうすると，相関のある場合，ない場合の両方に対応することが可能になります．

は直方体領域で与えられるのではなく，Ｆ分布によって規定される楕円体領域となることを示しました．ちなみに，パラメータの信頼区間が正規分布ならば，対応する母数空間はχ^2分布で規定される楕円体領域となります．直方体は楕円体領域の近似解にすぎないのです．

　そして，相関係数が０に近いときには長方形が第０近似解となるような太った楕円が与えられますが，相関が大きいほど細長い楕円になりますから，相関が考慮されることによって，誤差は一回り減少することがわかります．さて，この直観的な幾何学的表現を，もう少し解析的に評価してみることにしましょう．

　相関を考慮しない場合の信頼楕円は，半径√σiiの楕円体に相当しますから，その体積は分散・共分散行列Δの主対角要素の積の平方根のｍ乗

　一方，分散・共分散行列の固有値をλ1，・・・，λmとすると，楕円体の半径は√λiになりますから，その体積は

　分散・共分散行列は正定値行列ですから，すべての固有値は正であり，また，分散・共分散行列の行列式の値は

　行列式の値|Δ|は，ｍ＝２なら平行四辺形の面積，ｍ＝３なら平行六面体の体積となります．これはｍ次元単体の体積のｍ！倍です．もし，推論

が本当なら，教科書に載っているはずです．そこで，手持ちの本をあたっていたところ「線形作用素への誘い」古田孝之著（培風館）ｐ２０に，「アダマールの定理」を発見しました．

　それによると，上の推論は表現形は若干違うものの，アダマールの定理そのものです．アダマールの定理は，平行六面体の体積はノルムの積によって上から抑えられるという非常に興味深い事実を示しているのであって，ベクトルの内積はノルムの積より小さいというシュワルツの不等式の拡張であると見なすこともできます．

【参】すべての辺の長さが等しい平行六面体格子（菱形体格子）をつくってみると，辺が互いの６０°の角度をなすようにすると，平行六面体の体積は最小値となる．

　次に，標本空間において，最小２乗モデル式：ｙ＝ｆ（ｘ，θ）の推定誤差Δｙについて考えてみることにしましょう．

　この際，パラメータ数が３の場合についてのみ述べますが，「誤差伝播の法則」（propagation of error）を使って，

　この式は（本質的なところだけを抜き出すと）ピタゴラスの定理にほかならないので，母数空間が直方体領域で与えられるという誤解を招きますし，また，実際例に応用してみると誤差が大きすぎて，実用性には乏しいことがわかりました．なぜそうなるかというと，それは母数間の相関を考えていないからです．

　そこで，もう少しつじつまがあっていて，しかもいろいろな場面の応用できるましな方法に改良する余地ありということで，ピタゴラスの定理を一般化して得られる余弦定理にあたる誤差公式を考えてみることが必要になりました．

　最近の誤差解析の本には，やっとこの形の式が取り上げられるようになってきましたが，この改良式は母数間の相関を考慮したもので，余弦定理というよりも楕円体そのものを表す式となっています．すなわち，２つ（３つ）の母数がある場合，その同時信頼区間は長方形（直方体）領域で与えられるものではなく，楕円（楕円体）となることを示しているのであって，幾何学的には２次元ならば長方形に内接する楕円，３次元ならば直方体に内接する楕円体すなわちラグビーボール状態を考えることに相当します．

　　（Δｙ）^2＝（Δａ）^2＋（Δｂ）^2＋（Δｃ）^2＋２ΔａΔｂ＋２ΔｂΔｃ＋２ΔｃΔａ・・・・・（１）

では，前者のほうが小さいことが以下のようにして示されます（数学的な証明にはほど遠いのですが・・・）．

　厳密には「相関を考えたｔ乱数を発生させる」必要があっても，それは簡単にはいかないことがわかりました．可能なことことは，

ですが，通常，モンテカルロ法というと相関は考慮していないわけですから，相関を考えたｔ乱数を発生させる必要がある場合は，相関を考えた正規乱数を発生させることで代用すべきと思われます．1）は2），3）よりも優れた結果を与えてくれます．

　大抵のプログラム言語は乱数発生関数をもっていますが，ｔ分布など特定の分布に従うものはありません．一方，Mathematicaの標準パッケージでは分布関数とそれに従う乱数を発生させることができ，StudentDistributionがｔ分布にに相当します．しかし，多変量ｔ分布のように，標準パッケージ組み込み以外の分布関数について乱数を発生させことはできません．いまのところ，多変量正規乱数，多変量ｔ乱数などに対しては乱数発生アルゴリズムを自ら組み込むしかなさそうです．

　一方，誤差を見積もる問題は，シミュレーションに頼らなくとも，誤差伝播の法則を使って解決できることがわかりました．誤差伝播の法則は，分布形に関わらず（厳密にはコーシー分布を除いて）成り立ちますし，相関を考慮することも簡単です．相関を考慮すると，考慮しない場合より誤差が小さくなることは本文中で証明したとおりです．また，以前に，積分値の誤差は誤差の積分値よりも小さくなることを証明しましたが，これらの点でも，誤差伝播の法則は利用価値の高いものになっています．

　特性関数は積率母関数同様に和の分布を求めるときなど利用されていますが，ここでは，区間（０，１）の一様分布の特性関数が

となることを利用して，一様乱数ｒi（０～１）をｎ個合計したものの分布が，ｎ→∞の極限で正規分布になることを示してみましょう．

の両辺を比較することにより，Σ1/k2=π2/6，Σ1/k4=π4/90が計算される．Σ1/k2はゼータ関数ζ(2)に，Σ1/k4はゼータ関数ζ(4)に相当する．

　正規分布Ｎ（μ，σ^2）の特性関数はexp(iμt-σ2t2/2)ですから，この結果はｎ個の独立した一様乱数の和の分布は平均値ｎ／２，分散ｎ／１２の正規分布に近づくことを示しています．

　次に，一様乱数をもとに正規乱数を発生させる具体的な方法を述べてみましょう．中心極限定理によって，区間（０，１）の一様乱数ｒi（平均１／２，分散１／１２）をｎ個合計したものの分布は平均値μがｎ／２，分散σ^2がｎ／１２の正規分布に近くなりますから，正規化して

実際，一様分布に従う確率変数１２個ずつの平均をとり１００個のデータから構成したヒストグラムは元の一様分布とは似ても似つかない滑らかな分布となります．

　なお，中心極限定理を利用した正規乱数発生法は，１２個の一様乱数から１個の正規乱数が得られる効率のよくない方法で，それに比べ，ボックス・ミューラー法はかなり効率がよくなっています．

　z1,z2が正規分布N(0,1)にしたがい，独立のとき(z1^2+z2^2)^(1/2)はレイリー分布にしたがいます．したがって，レイリー乱数を発生させることができると正規乱数に変換することができます．ｒを一様乱数とすると(-2lnr)は平均値２の指数分布，したがって(-2lnr)^(1/2)はレイリー乱数，また2πrは(0,2π)の一様乱数になります．

　r^2=-2logr1,θ=2πr2ですから，これにより，２個の一様乱数ｒ1，ｒ2から互いに独立に標準正規分布に従う２個の正規乱数ｚ1，ｚ2

を作りだすことができます．このように，ボックス・ミューラー法は，標的問題の解であるレイリー分布を応用していると考えることができます．