以下はRmdファイルです。カイ二乗分布の多次元拡張版であることこの2点を中心にWishart分布について確認する異にする。あり、分布の形状を決めている。と1対1対応できるから、Wishart分布は多変量正規分布の形状決定因子に関する確率密度分布という位置づけになる。カイ二乗分布の多次元拡張版であることカイ二乗分布は標準正規分布の変数の値の二乗に関する確率密度分布である。カイ二乗分布に従う確率変数を1行1列の行列であるとみなし、これを正方行列に一般化する。確率変数として扱える。確率変数に確率密度を持たせると、サポート全体に確率密度関数が取れるが、それがWishart分布である。このようなWishart分布は、1行1列のときに自由度1のカイ二乗分布が一変数・一次元の正規分布に対応したように、何かしらの多変量正規分布と対応づく。サポートとするという点で多変量正規分布の形状決定変数である分散共分散行列と結びついていたように、この対応関係の存在も素直に納得できるものかもしれない。この2点に注目しつつ、Wishart分布をいじってみる、というのが、この文書の目的である。Rの準備その前に必要なパッケージなどの準備をしておく。サンプリングをし、そのサンプルがどのような分布を作るかを考える。正規分布の中心は原点0に固定し、分布の形状にのみ着目することとする。k$行列とその行列から生成される何かしらが、どのような分布を取るかを考えるために、Nサンプリングセットを用意することとする。k$個の値をアレイXに納めることとする。以下の記載では、n,k,Nをこの意味合いで用いる。標準偏差```まずは標準正規分布に関して、この条件で実行してみる。1$行列の$(1,1)$成分は平均0、標準偏差Sの正規分布を作る。正規分布を用いて発生したランダムサンプルからカイ二乗分布に従う変数のサンプルに変換する。この変換を、この文書の中では、カイ二乗化・Wishart化とでも呼ぶことにする。N個のサンプリングセットのそれぞれについてカイ二乗化・Wishart化してN個のアイテムを作成する。関数で書けば次のようになる。f(X[i,,])$$今の場合は、値の二乗に変換するので、次のようにすればよい。X^2mean(Y[,1,1])hist(Y[,1,1])```このYが作る確率密度分布は自由度1のカイ二乗分布である。1次元標準正規分布に従う確率変数Xの二乗であるYは自由度1のカイ二乗分布に従う念のため、確かめておく。n=1,k=1に限らず一般化する。ここで、カイ二乗化は単に二乗するだけであった。このままでは、Wishart分布が行列の分布であるので不都合である。n=1,p=1の場合は、発生させたX[i,,]が今は長さ1のベクトルであり、それをカイ二乗化した後のY[i,,]も長さ1のベクトルであったものを、行列的に扱いたい、と言うことである。X[i,,]$$Rの関数として作成しておく。Z)}```と定義することにする。k$行列の分布について考えることができる。行列の要素を使って書いておく。自由度dのカイ二乗分布とはn=1,k=1のときに自由度1のカイ二乗分布が得られた。カイ二乗分布の自由度1はn=1とk=1のどちらから来ているかを考えることにする。k$行列である確率変数を生じるから、k=1であろう。いうことは自由度dのカイ二乗分布に従う確率変数をもたらすのは、n=d,k=1の場合だろうとあたりがつく。d次元標準正規分布はd個の値のそれぞれを1次元標準正規分布から発生させて、長さdのベクトルにして得ることができる。n=d(=3たとえば),k=1,N=1000でやってみよう。アレイXのすべての要素が標準正規乱数なので以下のように作ることができる。array(rnorm(N*n*k,m,S),c(N,n,k))```N個のサンプリングセットについてカイ二乗化・Wishart化計算を行うと以下のようになる。my.wishart(X[i,,])}```関数にしておこう。Y}```関数を用いてカイ二乗化・Wishart化し直しておく。my.wishart.set(X)dim(Y)hist(Y[,1,1])```これが確かに自由度k=dのカイ二乗分布に従っているかどうかを確認してみよう。rchisq(N,d)plot(sort(Y[,1,1]),sort(Y.))abline(0,1,col=2)```自由度dを変えて同じことをやってみる。元の正規分布の形を変える〜分散を変える〜これまではn=1の場合とn=dの場合とについて、標準偏差S=1の標準正規分布に従う確率変数Xに基づいて、そのカイ二乗化・Wishart化した確率変数Yの分布が自由度nのカイ二乗分布に従うことを見てきた。Sの値を変えるとどうなるのかを確かめてみる。rchisq(N,d)plot(sort(Y[,1,1]),sort(Y.))abline(0,1,col=2)```このように、自由度dのカイ二乗分布には乗らない。次のようにすれば、Yのあ体は$S^2$倍になっていることがわかる。```{r}plot(sort(Y[,1,1])/(S^2),sort(Y.))abline(0,1,col=2)```念のため、dの値を変えて同じことをやってみる。分散が不明な正規分布からのデータを得たときSの値が変わると、カイ二乗化・Wishart化した値として得られやすい値が変わることがわかった。このことを利用して、標準偏差Sの正規分布からのサンプルデータであることがわかっている値がn=d個のあったとき、そのカイ二乗か・Wishart化した値を計算すれば、元の正規分布のSの値がいくつであったのかを尤度によって推定することができるだろう。my.wishart.set(X)Y```この1つのYの値だけから考えようというわけである。Sの値の候補をたくさん挙げて、そのそれぞれについて、Yの値が得られる確率を計算すればそれが尤度である。Yの値を$\frac{1}{S^2}$すれば自由度nのカイ二乗分布に従っているはずだから、自由度nのカイ二乗分布の確率密度を計算すればよいだろう。S.candidates[which(Like==max(Like))]max.S.candidateSplot(S.candidates,Like,type="l")abline(v=max.S.candidate,col=2)abline(v=S,col=3)```データセットが1つだけではあるが、n=10個の値があるのでそこそこ良い推定値になっているといったところだろう。これがd=1でやると、Yの値が正の値である限り、Sの値が大きければ大きいほど尤度が上がるので、この方法での推定はうまく行かない。それを示したのが以下である。多変量正規分布に拡張するこれまでは、一度にサンプリングする数nが1であったり、そうでなかったりするものの、サンプリングのもととなる分布は正規分布であった。サンプルから元の正規分布のSを尤度に照らして推定することをやってきた。今度は、サンプリングのもととなる分布を多変量正規分布に拡張する。以下のようにrnorm()関数をrmvnorm()関数に換えること、それに伴い、Sをスカラー値から分散共分散行列Sigmaに換えることでデータ作成は容易にできる。多変量正規分布の分散共分散行列と楕球と回転カイ二乗化・Wishart化に話を進める前に、多変量正規分布とSigmaとの関係について確認しておくことにする。Sigmaは分散共分散行列であるので、どんな行列でもよいというわけではない。V$と書き表せる必要がある。Vは回転行列であるので$V^{-1}=V^T$であることにも注意する。このような性質を持つ行列を半正定値行列と呼ぶ。うち最も標準的なそれは、変量間に相関がなく、すべての変量の分散が1であるようなそれである。そのような分布に等高線(等値線)を引くと多次元球面になる。変量の数をk=2としてやってみる。対応する$\Sigma$は単位行列(対角成分がすべて1で非対角成分がすべて0)である。今、すべての変量間に相関はないが、各々の変量の分散が1とは限らないときの$\Sigma$は、スケール行列と呼ばれるものとなる。スケール行列は、その対角成分が任意の非負の値であって、非対角成分は0であるような行列である。このような分布に等値線を引くと楕球になり、その軸は各変量に対応した軸となる。V$と分解すると、Uはこのスケール行列そのものであり、Vはもう勝手に回転させてはいけないので単位行列(もしくは反転に相当するような回転を表した行列)となる。固有値分解eigen(Sigma)lambdas```今、Uをスケール行列とし、Vを回転行列にすると、それに対応する多変量正規分布の等値線は、楕球を回転したものとなる(相変わらず楕球であるが、その軸が変量の軸と一致していないものとなる)。分散共分散行列についての確認が終わったので、カイ二乗化・Wishart化してみることにする。カイ二乗分布のときに見たように、あるサンプルデータから、そのサンプルが取られた元の分布形状に興味があるとする。そのようなとき(多変量正規分布の形を知りたいとき)、そこそこの数のデータが無ければ推定ができないことはわかる。k$くらいはなくてはならない。そのようなnを設定しよう。k=2の場合で、2変量が相互に独立で分散がどちらも1であるような場合から始める。2$行列であるカイ二乗化・Wishart化のYの(1,1),(2,2)成分は100付近でばらつき、(1,2),(2,1)成分は0付近でばらついていることがわかる。nの値を変えてみよう。my.wishart.set(X)apply(Y,c(2,3),mean)```カイ二乗化・Wishart化のYの(1,1),(2,2)成分はおよそ30、(1,2),(2,1)成分はおよそ0であることがわかる。(1,1),(2,2)成分はnのようだ。lambdas```では、$\Sigma$の対角成分を1からずらしてみることにする。Yの(1,1),(2,2)成分は対角成分(2,0.3)にnを掛けた値になっているようだ。念のためnと対角成分の値を変えてやってみる。lambdas```ここまでは$\Sigma$の非対角成分が0であった。回転成分Vを加えて変量間に関連を持たせてみる。lambdas```今度はYの(1,2),(2,1)成分が0でなくなり、(1,1),(2,2)成分も固有値のn倍とはなっていない。Sigmaapply(Y,c(2,3),mean)```代わりに、$\Sigma$のn倍になっているようだ。念のためn,kを変えてもう一度やっておく。\Sigma$を平均として、その周辺にばらついていることがわかった。いうことはYを知ることで$\Sigma$について推定ができることになる。1$のときに、多変量正規分布の$\Sigma$をサンプル値からYを計算し、それを何かに照らして$\Sigma$に対する尤度を計算することで、$\Sigma$の推定をしたい。k$行列であり、その成分が$k^2$個あるので、それを網羅するのはちょっと面倒である。k=2に限定すると、2個の固有値と回転角との3変数で分散共分散行列を指定することができるから、この3つを変数に網羅的に$\Sigma$の候補を作成し、それに対するWishart分布の確率を計算してみる。k=2の場合において、適当に$\Sigma$を作成し、それに基づく観測データセットを作り、Yを計算することとする。my.wishart.set(X)Y[1,,]```ついで、その上で、そのようなYを得る確率をさまざまな$\Sigma$候補行列に関してWishart分布の確率密度を計算してみる。$\Sigma$候補は次のように、2つの固有値lambda1,2と回転角thetaとの3変数で網羅的に作成することとする。lambda1,2とthetaから分散共分散行列を作る関数を作成する。V}```今、lambda1は$\Sigma$の第1の固有値に固定し、lambda2とthetaを変化させてみる。params[,1]*sin(params[,2])plot(x.coord,y.coord,col=rgb(1-st.Like,st.Like,1),pch=20)```次にlambda2を固定してみる。params[,1]*sin(params[,2])plot(x.coord,y.coord,col=rgb(1-st.Like,st.Like,1),pch=20)```最後にthetaを固定してみる。多変量正規分布の分散共分散行列をランダムに発生させるとき〜BUGS〜多変量正規分布における分散共分散行列$\Sigma$の幾何的な意味は確認した。\Sigma$を中心とした確率密度分布を持つことも確認した。従って多変量正規分布からのサンプルセットがあったときに、元の多変量正規分布の$\Sigma$を尤度に基づいて推定する、という考え方があることも確認した。今度は、半正定値行列を色々に発生させたいときに、Wishart分布からのランダムな変数(行列)を発生させる、という話をする。半正定値行列をランダムに発生させることは、多変量正規分布の分散共分散行列を推定したいときに、BUGSを回す時に必要となる。RのMCMCpackにあるrwish()関数は、自由度nと半正定値行列Zを引数としてとり、半正定値行列をランダム変数(行列)として返す関数である。\Sigma$であるような変数(行列)が生じる。