BUGSのときのWishart分布

• ここのところBUGSをやっていて、こちらにWishart分布についてメモした
• その続き
• BUGSで多変量正規分布の分散共分散行列を推定したいとき、なるべくばらばらな半正定値行列をランダムに発生させたい
• そのときにWishart分布から半正定値行列を発生させるわけだけれど、Wishart分布からのランダム行列発生には(1)自由度と(2)発生させる半正定値行列の期待値行列の2つの引数を与える必要があるのだが
• それは、多変量間に相関があるのかないのか不明なとき(ないかもしれないし正の相関があるかもしれないし負のそれがあるかもしれないとき)には、対角行列であって、各変量の分散として適当な値を選び(何かは選ばなくてはならな、本当に何の情報もなければ1にするのもよい)の上で、ばらつきを大きくするために自由度は変量の数にするのがよい、とBUGSのtutorial等に書いてある
• その通りなのだけれど、それってどういうこと？というのを理解するために書いたのが以下のRmarkdonw文書
• 以下はRmdファイルです。こちらなどを参考にすれば、自由にhtml文書化・epub化できます(が、それが面倒な場合には、kindleで)
• こちらから

Wishart分布?BUGSのために?

# Wishart 分布

Wishart分布には、大きく２つの特徴がある
* 対称性半正定値行列全体に定義された確率密度分布であること
* カイ二乗分布の多次元拡張版であること

この２点を中心にWishart分布について確認する異にする。

## 対称性半正定値行列全体をサポートすること

多変量正規分布の分散共分散行列は対称性半正定値行列であり、分布の形状を決めている。
したがって、中心を固定すれば、多変量正規分布全体は対称性半正定値行列と１対１対応できるから、Wishart分布は多変量正規分布の形状決定因子に関する確率密度分布という位置づけになる。

## カイ二乗分布の多次元拡張版であること

カイ二乗分布は標準正規分布の変数の値の二乗に関する確率密度分布である。
カイ二乗分布に従う確率変数を１行１列の行列であるとみなし、これを正方行列に一般化する。その一般化にあたって、対称性半正定値正方行列に拡張すると、その対称性半正定値行列が確率変数として扱える。
確率変数に確率密度を持たせると、サポート全体に確率密度関数が取れるが、それがWishart分布である。

このようなWishart分布は、１行１列のときに自由度１のカイ二乗分布が一変数・一次元の正規分布に対応したように、何かしらの多変量正規分布と対応づく。
すでにWishart分布は対称性半正定値行列をサポートとするという点で多変量正規分布の形状決定変数である分散共分散行列と結びついていたように、この対応関係の存在も素直に納得できるものかもしれない。

この２点に注目しつつ、Wishart分布をいじってみる、というのが、この文書の目的である。

## Rの準備

その前に必要なパッケージなどの準備をしておく。

```{r}
library(mvtnorm) # 多変量正規分布を扱うため
library(MCMCpack) # Wishart分布を扱うため
library(GPArotation) # 多次元回転行列を扱うため
library(rgl) # ３次元プロットのため
```

## 分布を取るための準備

今から正規分布・多変量正規分布からサンプリングをし、そのサンプルがどのような分布を作るかを考える。

また、正規分布の中心は原点0に固定し、分布の形状にのみ着目することとする。

そのために、k-変量正規分布から、n個の相互に独立なサンプルをとり、それが作る$n\times k$行列を1サンプリングセットとし、そのような$n\times k$行列とその行列から生成される何かしらが、どのような分布を取るかを考えるために、Nサンプリングセットを用意することとする。

この$N\times n \times k$個の値をアレイXに納めることとする。

以下の記載では、n,k,Nをこの意味合いで用いる。

## いわゆる標準正規分布と自由度１のカイ二乗分布

```{r}
n <- 1 # １データセットを構成する
k <- 1 # 空間の次元
N <- 1000 # 分布を取るために発生させるランダムサンプルの数
m <- rep(0,k) # 原点座標
S <- 1 # 標準偏差
```

まずは標準正規分布に関して、この条件で実行してみる。

```{r}
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rnorm(n,m,S)
}
hist(X[,1,1]) # 正規分布
```

このようにして発生したN個の$1\times 1$行列の$(1,1)$成分は平均0、標準偏差Sの正規分布を作る。

```{r}
mean(X[,1,1])
sd(X[,1,1])
```

## カイ二乗化・Wishart化

### n=1,k=1の場合
正規分布を用いて発生したランダムサンプルからカイ二乗分布に従う変数のサンプルに変換する。

この変換を、この文書の中では、カイ二乗化・Wishart化とでも呼ぶことにする。

N個のサンプリングセットのそれぞれについてカイ二乗化・Wishart化してN個のアイテムを作成する。

関数で書けば次のようになる。

$$
Y_i = f(X[i,,])
$$

今の場合は、値の二乗に変換するので、次のようにすればよい。
```{r}
Y <- X^2
mean(Y[,1,1])
hist(Y[,1,1])
```

このYが作る確率密度分布は自由度1のカイ二乗分布である。

* 1次元標準正規分布に従う確率変数Xの二乗であるYは自由度1のカイ二乗分布に従う
念のため、確かめておく。
```{r}
Y. <- rchisq(N,1)
plot(sort(Y),sort(Y.))
abline(0,1,col=2)
```

### n=1,k=1に限らず一般化する。

ここで、カイ二乗化は単に二乗するだけであった。
このままでは、Wishart分布が行列の分布であるので不都合である。

n=1,p=1の場合は、発生させたX[i,,]が今は長さ1のベクトルであり、それをカイ二乗化した後のY[i,,]も長さ1のベクトルであったものを、行列的に扱いたい、と言うことである。

長さ1のベクトルは$1 \times 1$行列とみなせるから、そのようにするとして、カイ二乗化を$n \times k$行列$X[i,,]$を用いて

$$
Y_i = f(X[i,,]) = (X[i,,])^T X[i,,]
$$
Rの関数として作成しておく。

```{r}
my.wishart <- function(Z){
  as.matrix(t(Z) %*% Z)
}
```
と定義することにする。
これにより$k \times k$行列がN個得られるから、$k \times k$行列の分布について考えることができる。

行列の要素を使って書いておく。

$$
\begin{equation}
Y_i = (y^i_{u,v}); u,v = 1,2,...,k\\
y^i_{u,v} = \sum_{z=1}^n X[i,z,u]X[i,z,v]
\end{equation}
$$


## 自由度dのカイ二乗分布とは

n=1,k=1のときに自由度1のカイ二乗分布が得られた。

カイ二乗分布の自由度1はn=1とk=1のどちらから来ているかを考えることにする。

自由度dのカイ二乗分布に従う確率変数は、スカラー($1 \times 1$ 行列)な変数である。
カイ二乗か・Wishart化では、$k\times k$行列である確率変数を生じるから、k=1であろう。
ということは自由度dのカイ二乗分布に従う確率変数をもたらすのは、n=d,k=1の場合だろうとあたりがつく。

やってみよう。

d次元標準正規分布はd個の値のそれぞれを1次元標準正規分布から発生させて、長さdのベクトルにして得ることができる。

n=d(=3たとえば),k=1,N=1000でやってみよう。

アレイXのすべての要素が標準正規乱数なので以下のように作ることができる。


```{r}
d <- 3
n <- d
k <- 1
N <- 1000
m <- rep(0,k) # 原点座標
S <- 1 # 標準偏差
X <- array(rnorm(N*n*k,m,S),c(N,n,k))
```

N個のサンプリングセットについてカイ二乗化・Wishart化計算を行うと以下のようになる。
```{r}
Y <- array(0,c(N,k,k))
for(i in 1:N){
  Y[i,,] <- my.wishart(X[i,,])
}
```

関数にしておこう。

```{r}
my.wishart.set <- function(X){
  d <- dim(X)
  Y <- array(0,c(d[1],d[3],d[3]))
  for(i in 1:d[1]){
    Y[i,,] <- my.wishart(X[i,,])
  }
  Y
}
```

関数を用いてカイ二乗化・Wishart化し直しておく。
```{r}
Y <- my.wishart.set(X)
dim(Y)
hist(Y[,1,1])
```

これが確かに自由度k=dのカイ二乗分布に従っているかどうかを確認してみよう。

```{r}
Y. <- rchisq(N,d)
plot(sort(Y[,1,1]),sort(Y.))
abline(0,1,col=2)
```

自由度dを変えて同じことをやってみる。

```{r}
d <- 8
n <- d
k <- 1
N <- 1000
m <- rep(0,k) # 原点座標
S <- 1 # 標準偏差
X <- array(rnorm(N*n*k,m,S),c(N,n,k))
Y <- my.wishart.set(X)
hist(Y)
Y. <- rchisq(N,d)
plot(sort(Y[,1,1]),sort(Y.))
abline(0,1,col=2)
```

## 元の正規分布の形を変える〜分散を変える〜

これまではn=1の場合とn=dの場合とについて、標準偏差S=1の標準正規分布に従う確率変数Xに基づいて、そのカイ二乗化・Wishart化した確率変数Yの分布が自由度nのカイ二乗分布に従うことを見てきた。

Sの値を変えるとどうなるのかを確かめてみる。

```{r}
d <- 1
n <- d
k <- 1
N <- 1000
m <- rep(0,k) # 原点座標
S <- 3 # 標準偏差
X <- array(rnorm(N*n*k,m,S),c(N,n,k))
Y <- my.wishart.set(X)
hist(Y)
Y. <- 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の値を変えて同じことをやってみる。
```{r}
d <- 5
n <- d
k <- 1
N <- 1000
m <- rep(0,k) # 原点座標
S <- 3 # 標準偏差
X <- array(rnorm(N*n*k,m,S),c(N,n,k))
Y <- my.wishart.set(X)
hist(Y)
Y. <- rchisq(N,d)
par(mfcol=c(1,2))
plot(sort(Y[,1,1]),sort(Y.))
abline(0,1,col=2)
plot(sort(Y[,1,1])/(S^2),sort(Y.))
abline(0,1,col=2)
par(mfcol=c(1,1))
```

## 分散が不明な正規分布からのデータを得たとき

Sの値が変わると、カイ二乗化・Wishart化した値として得られやすい値が変わることがわかった。

このことを利用して、標準偏差Sの正規分布からのサンプルデータであることがわかっている値がn=d個のあったとき、そのカイ二乗か・Wishart化した値を計算すれば、元の正規分布のSの値がいくつであったのかを尤度によって推定することができるだろう。

やってみる。

```{r}
S <- runif(1)*10 # ランダムに決める
d <- 10
n <- d
k <- 1
N <- 1 # データセットは１個だけ
m <- rep(0,k) # 原点座標

X <- array(rnorm(N*n*k,m,S),c(N,n,k))
Y <- my.wishart.set(X)
Y
```

この1つのYの値だけから考えようというわけである。

Sの値の候補をたくさん挙げて、そのそれぞれについて、Yの値が得られる確率を計算すればそれが尤度である。
Yの値を$\frac{1}{S^2}$すれば自由度nのカイ二乗分布に従っているはずだから、自由度nのカイ二乗分布の確率密度を計算すればよいだろう。

```{r}
S.candidates <- seq(from=0,to=20,length=100)
S.candidates <- S.candidates[-1] # S=0はあり得ないので除く
Like <- dchisq(Y/(S.candidates^2),n)
max.S.candidate <- S.candidates[which(Like==max(Like))]
max.S.candidate
S
plot(S.candidates,Like,type="l")
abline(v=max.S.candidate,col=2)
abline(v=S,col=3)
```

データセットが1つだけではあるが、n=10個の値があるのでそこそこ良い推定値になっているといったところだろう。

これがd=1でやると、Yの値が正の値である限り、Sの値が大きければ大きいほど尤度が上がるので、この方法での推定はうまく行かない。
それを示したのが以下である。

```{r}
S <- runif(1)*10 # ランダムに決める
d <- 1
n <- d
k <- 1
N <- 1 # データセットは１個だけ
m <- rep(0,k) # 原点座標

X <- array(rnorm(N*n*k,m,S),c(N,n,k))
Y <- my.wishart.set(X)
Y

S.candidates <- seq(from=0,to=20,length=100)
S.candidates <- S.candidates[-1] # S=0はあり得ないので除く
Like <- dchisq(Y/(S.candidates^2),n)
max.S.candidate <- S.candidates[which(Like==max(Like))]
max.S.candidate
S
plot(S.candidates,Like,type="l")
abline(v=max.S.candidate,col=2)
abline(v=S,col=3)
```

## 多変量正規分布に拡張する

これまでは、一度にサンプリングする数nが1であったり、そうでなかったりするものの、サンプリングのもととなる分布は正規分布であった。
そしてサンプルから元の正規分布のSを尤度に照らして推定することをやってきた。

今度は、サンプリングのもととなる分布を多変量正規分布に拡張する。
以下のようにrnorm()関数をrmvnorm()関数に換えること、それに伴い、Sをスカラー値から分散共分散行列Sigmaに換えることでデータ作成は容易にできる。

```{r}
n <- 1000
k <- 2
N <- 1
m <- rep(0,k) # 原点座標
Sigma <- diag(rep(1,k)) # 分散共分散行列。
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
plot(X[1,,],asp=1)
Y <- my.wishart.set(X)
Y[1,,]
```

## 多変量正規分布の分散共分散行列と楕球と回転

カイ二乗化・Wishart化に話を進める前に、多変量正規分布とSigmaとの関係について確認しておくことにする。

Sigmaは分散共分散行列であるので、どんな行列でもよいというわけではない。
対称行列であって、固有値分解すると、回転行列Vと非負の固有値を対角成分とする対角行列$U$とを用いて$\Sigma=V^T U V$と書き表せる必要がある。Vは回転行列であるので$V^{-1}=V^T$であることにも注意する。

このような性質を持つ行列を半正定値行列と呼ぶ。

ところで、このV,Uと多変量正規分布の関係を考えることとする。

多変量正規分布のうち最も標準的なそれは、変量間に相関がなく、すべての変量の分散が1であるようなそれである。

そのような分布に等高線(等値線)を引くと多次元球面になる。

変量の数をk=2としてやってみる。
```{r}
n <- 1000
k <- 2
N <- 1
m <- rep(0,k) # 原点座標
Sigma <- diag(rep(1,k)) # 分散共分散行列。
Sigma
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
plot(X[1,,],asp=1)
```

また、そのような標準多変量正規分布に対応する$\Sigma$は単位行列(対角成分がすべて1で非対角成分がすべて0)である。

今、すべての変量間に相関はないが、各々の変量の分散が1とは限らないときの$\Sigma$は、スケール行列と呼ばれるものとなる。スケール行列は、その対角成分が任意の非負の値であって、非対角成分は0であるような行列である。

このような分布に等値線を引くと楕球になり、その軸は各変量に対応した軸となる。

このような$\Sigma$を$V^T U V$と分解すると、Uはこのスケール行列そのものであり、Vはもう勝手に回転させてはいけないので単位行列(もしくは反転に相当するような回転を表した行列)となる。

```{r}
n <- 1000
k <- 2
N <- 1
m <- rep(0,k) # 原点座標
lambdas <- c(2,0.5) # 対角成分
Sigma <- diag(lambdas) # 分散共分散行列。
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
plot(X[1,,],asp=1)
Sigma
# 固有値分解
eigen(Sigma)
lambdas
```

今、Uをスケール行列とし、Vを回転行列にすると、それに対応する多変量正規分布の等値線は、楕球を回転したものとなる(相変わらず楕球であるが、その軸が変量の軸と一致していないものとなる)。

```{r}
n <- 1000
k <- 2
N <- 1
m <- rep(0,k) # 原点座標
lambdas <- c(2,0.5) # 対角成分
# 回転行列
V <- Random.Start(k)
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
plot(X[1,,],asp=1)
Sigma
eigen(Sigma)
lambdas
```

## 多変量正規分布からのサンプルとカイ二乗化・Wishart化

多変量正規分布と分散共分散行列についての確認が終わったので、カイ二乗化・Wishart化してみることにする。

カイ二乗分布のときに見たように、あるサンプルデータから、そのサンプルが取られた元の分布形状に興味があるとする。

そのようなとき(多変量正規分布の形を知りたいとき)、そこそこの数のデータが無ければ推定ができないことはわかる。

$n \ge k$くらいはなくてはならない。そのようなnを設定しよう。

まずはk=2の場合で、2変量が相互に独立で分散がどちらも1であるような場合から始める。

```{r}
n <- 100
k <- 2
N <- 1000
m <- rep(0,k) # 原点座標
lambdas <- rep(1,k) # 対角成分
# 回転行列
# V <- Random.Start(k)
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
Y <- my.wishart.set(X)
apply(Y,c(2,3),mean)
boxplot(cbind(Y[,1,1],Y[,1,2],Y[,2,1],Y[,2,2]))
```

$$\Sigma = \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}$$のときにn=100で実施すると、$2\times 2$行列であるカイ二乗化・Wishart化のYの(1,1),(2,2)成分は100付近でばらつき、(1,2),(2,1)成分は0付近でばらついていることがわかる。

nの値を変えてみよう。
```{r}
n <- 30
k <- 2
N <- 1000
m <- rep(0,k) # 原点座標
lambdas <- rep(1,k) # 対角成分
# 回転行列
V <- diag(rep(1,k))
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
Y <- 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のようだ。
```{r}
n <- 100
k <- 2
N <- 1000
m <- rep(0,k) # 原点座標
lambdas <- c(2,0.3) # 対角成分
# 回転行列
V <- diag(rep(1,k))
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
Y <- my.wishart.set(X)
apply(Y,c(2,3),mean)
n * lambdas
```
では、$\Sigma$の対角成分を1からずらしてみることにする。

Yの(1,1),(2,2)成分は対角成分(2,0.3)にnを掛けた値になっているようだ。
念のためnと対角成分の値を変えてやってみる。
```{r}
n <- 40
k <- 2
N <- 1000
m <- rep(0,k) # 原点座標
lambdas <- c(1.5,9) # 対角成分
# 回転行列
V <- diag(rep(1,k))
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
Y <- my.wishart.set(X)
apply(Y,c(2,3),mean)
n * lambdas
```

ここまでは$\Sigma$の非対角成分が0であった。
回転成分Vを加えて変量間に関連を持たせてみる。
```{r}
n <- 100
k <- 2
N <- 1000
m <- rep(0,k) # 原点座標
lambdas <- c(1.5,9) # 対角成分
# 回転行列
V <- Random.Start(k)
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
Sigma
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
Y <- my.wishart.set(X)
apply(Y,c(2,3),mean)
boxplot(cbind(Y[,1,1],Y[,1,2],Y[,2,1],Y[,2,2]))
n * lambdas
```

今度はYの(1,2),(2,1)成分が0でなくなり、(1,1),(2,2)成分も固有値のn倍とはなっていない。


```{r}
n * Sigma
apply(Y,c(2,3),mean)
```
代わりに、$\Sigma$のn倍になっているようだ。
念のためn,kを変えてもう一度やっておく。

```{r}
n <- 30
k <- 4
N <- 1000
m <- rep(0,k) # 原点座標
lambdas <- runif(k)*3 # 対角成分
# 回転行列
V <- Random.Start(k)
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
Sigma
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
Y <- my.wishart.set(X)
apply(Y,c(2,3),mean)
n * Sigma
```

結局、カイ二乗化・Wishart化のYは$n \times \Sigma$を平均として、その周辺にばらついていることがわかった。

ということはYを知ることで$\Sigma$について推定ができることになる。


## $\Sigma$が不明な多変量正規分布からのデータを得たとき

k=1のときに、正規分布のSをサンプル値からYを計算し、それをカイ二乗分布に照らしてSの値に対する尤度を計算することでSの推定をしたのと同じことを、$k\ge 1$のときに、多変量正規分布の$\Sigma$をサンプル値からYを計算し、それを何かに照らして$\Sigma$に対する尤度を計算することで、$\Sigma$の推定をしたい。

このときk=1の場合のカイ二乗分布に対応するのが$k \times k$の行列(対称半正定値行列)をサポートとする確率密度分布であるWishart分布である。

分散共分散行列は$k\times k$行列であり、その成分が$k^2$個あるので、それを網羅するのはちょっと面倒である。

k=2に限定すると、2個の固有値と回転角との3変数で分散共分散行列を指定することができるから、この3つを変数に網羅的に$\Sigma$の候補を作成し、それに対するWishart分布の確率を計算してみる。

まず、k=2の場合において、適当に$\Sigma$を作成し、それに基づく観測データセットを作り、Yを計算することとする。


```{r}
n <- 100
k <- 2
N <- 1
m <- rep(0,k) # 原点座標
lambdas <- runif(k)*2 # 対角成分
# 回転行列
the <- runif(1)*pi
V <- matrix(c(cos(the),sin(the),-sin(the),cos(the)),2,2)
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
Sigma
X <- array(0,c(N,n,k))
for(i in 1:N){
  X[i,,] <- rmvnorm(n=n, mean=m, sigma=Sigma)
}
Y <- my.wishart.set(X)
Y[1,,]
```

ついで、
その上で、そのようなYを得る確率をさまざまな$\Sigma$候補行列に関してWishart分布の確率密度を計算してみる。
$\Sigma$候補は次のように、2つの固有値lambda1,2と回転角thetaとの3変数で網羅的に作成することとする。

まず、lambda1,2とthetaから分散共分散行列を作る関数を作成する。
```{r}
# 3変数から2x2行列を作る関数
my.covM.2 <- function(lambda1,lambda2,theta){
  V <- matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)
  U <- diag(c(lambda1,lambda2))
  t(V) %*% U %*% V
}
```

今、lambda1は$\Sigma$の第1の固有値に固定し、lambda2とthetaを変化させてみる。
```{r}
lambda1 <- lambdas[1]
lambda2 <- seq(from=lambdas[2]/2,to=lambdas[2]*2,length=50)
theta <- seq(from=0,to=pi,length=50)
params <- expand.grid(lambda2,theta)
Like <- rep(0,length(params[,1]))
for(i in 1:length(params[,1])){
  S.candidate <- my.covM.2(lambda1,params[i,1],params[i,2])
  Like[i] <- dwish(Y[1,,],n,S.candidate) 
}
st.Like <- (Like-min(Like))/(max(Like)-min(Like))
x.coord <- params[,1]*cos(params[,2])
y.coord <- params[,1]*sin(params[,2])
plot(x.coord,y.coord,col=rgb(1-st.Like,st.Like,1),pch=20)
```
次にlambda2を固定してみる。
```{r}
lambda2 <- lambdas[2]
lambda1 <- seq(from=lambdas[1]/2,to=lambdas[1]*2,length=50)
theta <- seq(from=0,to=pi,length=50)
params <- expand.grid(lambda1,theta)
Like <- rep(0,length(params[,1]))
for(i in 1:length(params[,1])){
  S.candidate <- my.covM.2(params[i,1],lambda2,params[i,2])
  Like[i] <- dwish(Y[1,,],n,S.candidate) 
}
st.Like <- (Like-min(Like))/(max(Like)-min(Like))
x.coord <- params[,1]*cos(params[,2])
y.coord <- params[,1]*sin(params[,2])
plot(x.coord,y.coord,col=rgb(1-st.Like,st.Like,1),pch=20)
```

最後にthetaを固定してみる。
```{r}
theta <- the
lambda1 <- seq(from=lambdas[1]/2,to=lambdas[1]*2,length=50)
lambda2 <- seq(from=lambdas[2]/2,to=lambdas[2]*2,length=50)
params <- expand.grid(lambda1,lambda2)
Like <- rep(0,length(params[,1]))
for(i in 1:length(params[,1])){
  S.candidate <- my.covM.2(params[i,1],params[i,2],theta)
  Like[i] <- dwish(Y[1,,],n,S.candidate) 
}
st.Like <- (Like-min(Like))/(max(Like)-min(Like))
x.coord <- params[,1]
y.coord <- params[,2]
plot(x.coord,y.coord,col=rgb(1-st.Like,st.Like,1),pch=20)
```

## 多変量正規分布の分散共分散行列をランダムに発生させるとき〜BUGS〜

多変量正規分布における分散共分散行列$\Sigma$の幾何的な意味は確認した。

また、多変量正規分布からのサンプルセットがあったときに、そのカイ二乗化・Wishart化の$k\times k$行列Y(Scatter 行列と言う。そろそろ用語を出してもよいだろう)という確率変数は、半正定値行列をサポートとし、$n \times \Sigma$を中心とした確率密度分布を持つことも確認した。

従って多変量正規分布からのサンプルセットがあったときに、元の多変量正規分布の$\Sigma$を尤度に基づいて推定する、という考え方があることも確認した。

今度は、半正定値行列を色々に発生させたいときに、Wishart分布からのランダムな変数(行列)を発生させる、という話をする。

半正定値行列をランダムに発生させることは、多変量正規分布の分散共分散行列を推定したいときに、BUGSを回す時に必要となる。

RのMCMCpackにあるrwish()関数は、自由度nと半正定値行列Zを引数としてとり、半正定値行列をランダム変数(行列)として返す関数である。

このようなランダム変数(行列)をたくさん発生させると、平均が$n\times \Sigma$であるような変数(行列)が生じる。
以下で確かめる。

```{r}
n <- 30
k <- 4
lambdas <- runif(k)*3 # 対角成分
# 回転行列
V <- Random.Start(k)
Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
N <- 10000
R <- array(0,c(N,k,k))
for(i in 1:N){
  R[i,,] <- rwish(n, Sigma)
}

(mm <- apply(R,c(2,3),mean))
Sigma
mm/(Sigma*n)
```
```


多くの場合、多変量の間に正の相関があるか負の相関があるかについて中立の立場を取りたいので、分散共分散行列の非対角成分は正の値も負の値も取らせつつ、平均すれば0になるような半正定値行列を発生させたい。

従って、このような場合には、Wishart分布からのランダム変数(行列)発生にあたっては、対角行列を与えることが適切であることがわかる。

では、nが発生する変数(行列)の成分のばらつきに与える影響を確認するために、すべての対角成分を等しくし、さらに$n\tims \lambda$の値を1に固定した上で、nの値を変えて、発生する行列の成分のばらつきを調べてみることにする。

```{r}
k <- 4
ns <- seq(from=k,to=100,length=10)
res <- matrix(0,length(ns),k^2)
lambdas.pre <- rep(1,k)
for(j in 1:length(ns)){
  n <- ns[j]
  lambdas <- lambdas.pre/n # 対角成分
  # 回転行列
  V <- diag(rep(1,k))
  Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
  N <- 1000
  R <- array(0,c(N,k,k))
  for(i in 1:N){
   R[i,,] <- rwish(n, Sigma)
  }

  res[j,] <- apply(R,c(2,3),var)
  
}
matplot(ns,res,type="l",xlab="df",ylab="var")
abline(v=ns[1])
abline(h=0)
```

ここにみるように、自由度nが大きくなると分散は小さくなり0に収束していく様子がわかる。

また、rwish()関数にk未満のnを与えるとえらが出るのだが、それが示すように、分散を最大にするためにはn=kとすることが適当であることがわかる。

では、対角成分を1に揃えず適当に変えた場合にも分散とnとの関係が成り立つかどうかを確認する。

```{r}
k <- 4
ns <- seq(from=k,to=100,length=10)
res <- matrix(0,length(ns),k^2)
lambdas.pre <- runif(k)*5 # 対角成分を揃えない
for(j in 1:length(ns)){
  n <- ns[j]
  lambdas <- lambdas.pre/n # 対角成分
  # 回転行列
  V <- diag(rep(1,k))
  Sigma <- t(V) %*% diag(lambdas) %*% V # 分散共分散行列。
  N <- 1000
  R <- array(0,c(N,k,k))
  for(i in 1:N){
   R[i,,] <- rwish(n, Sigma)
  }

  res[j,] <- apply(R,c(2,3),var)
  
}
matplot(ns,res,type="l",xlab="df",ylab="var")
abline(v=ns[1])
abline(h=0)
```

対角成分がばらばらでも自由度が小さい方がばらつきが大きいことが確認できた。

従って、変量間に相関があるかどうかわからない場合の多変量正規分布のための半正定値行列のランダム発生においては、
非対角成分は0とし、対角成分自体は各変量の分散の期待値に合わせて選びつつ、なるべくばらばらの値を取る(それにより変量間の相関の強弱のばらつきが生み出される)ようにするためには、自由度nはkに一致させることが適当であることがわかる。