2017-07-15

逆Wishart分布

ベイズ推定正規分布 Wishart分布逆Wishart分布 χ二乗分布 Bartlett分解

正規分布推定をベイズで行うと、逆Wishart分布と言うのが出てくる
このあたりを整理する
Wishart分布
- 多変量正規分布の分散共分散行列を「平均」としてばらついている正方行列の分布

library(MCMCpack)
library(mvtnorm)
# 多変量正規分布の分散共分散行列
Sigma <- matrix(c(3,2,2,3),2,2)
Sigma
# 変量の個数
p <- dim(Sigma)[1]
# サンプル数
n <- 10
# p次元多変量正規分布から
# n標本を取って
# n標本間の内積行列t(X)Xを作ると
# それはもとの分散共分散行列と同じ次元の正方行列の形をしているが
# その個々の成分の期待値は
# もとの分散共分散の成分のn倍で
# その個々の成分の分散は
# 元の分散共分散行列の成分を使って、n(v_{ij}^2 + v_{ii}v_{jj})となる
# 以下、確かめる

n.iter <- 1000 # nサンプルをサンプリングする試行の繰り返し回数
series.XtX <- matrix(0,n.iter,p^2) # 出来るt(X)Xの成分を格納する箱
for(i in 1:n.iter){
  X <- rmvnorm(n,sigma=Sigma)
  tmp <- t(X) %*% X
  series.XtX[i,] <- c(tmp)
}
apply(series.XtX,2,mean) # 各成分の平均
c(n * Sigma)
apply(series.XtX,2,var) # 各成分の分散
c(n * (Sigma^2 + matrix(diag(Sigma),ncol=1) %*% matrix(diag(Sigma),nrow=1)))

- χ二乗分布の一般化

# p = 1、Sigma = 1のとき
# n個のサンプルをサンプリングすると
# t(X)Xは、n個の１次元標準正規乱数の距離の二乗の和
# それは自由度nのカイ二乗分布に従う
# 自由度nのカイ二乗乱数とのcoplotをすることでそれを確かめる
Sigma <- matrix(1,1,1)
# 変量の個数
p <- dim(Sigma)[1]
# サンプル数
ns <- 1:9
par(mfrow=c(3,3))
n.iter <- 1000 # nサンプルをサンプリングする試行の繰り返し回数
for(ii in ns){
  n <- ii
  series.XtX <- matrix(0,n.iter,p^2) # 出来るt(X)Xの成分を格納する箱
  for(i in 1:n.iter){
    X <- rmvnorm(n,sigma=Sigma)
    tmp <- t(X) %*% X
    series.XtX[i,] <- c(tmp)
  }
  plot(sort(series.XtX),sort(rchisq(n.iter,df=ii)),pch=20)
  abline(0,1,col=2)
}

- MCMCpack::rwish()関数はBartlett分解を利用してランダムWishart行列を作成している
  - その背景にあるのはBartlett分解
  - ある分散共分散行列Sigmaに従う正規乱数から得られるWishart分布からのランダム行列は：
    - $L = chol(Sigma)$ (Sigmaのコレスキー分解)のLを使って $LAA^TL^T$ として作れる
    - ただしAは対角成分が、自由度がp,p-1,...,1とだんだん変化するのカイ二乗乱数で、対角線下の三角形成分が標準正規乱数であるのようなもの→Bartlett分解の記事。
    - この分解も「カイ二乗分布の拡張からのランダムな行列A」と、それを「特定の分散共分散行列的なパターンにする」というように読めて、カイ二乗分布の一般化であることを感じさせる。

> rwish
function (v, S) 
{
    if (!is.matrix(S)) 
        S <- matrix(S)
    if (nrow(S) != ncol(S)) {
        stop(message = "S not square in rwish().\n")
    }
    if (v < nrow(S)) {
        stop(message = "v is less than the dimension of S in rwish().\n")
    }
    p <- nrow(S)
    CC <- chol(S)
    Z <- matrix(0, p, p)
    diag(Z) <- sqrt(rchisq(p, v:(v - p + 1)))
    if (p > 1) {
        pseq <- 1:(p - 1)
        Z[rep(p * pseq, pseq) + unlist(lapply(pseq, seq))] <- rnorm(p * 
            (p - 1)/2)
    }
    return(crossprod(Z %*% CC))
}

逆Wishart分布
- Wishart分布に従う行列があったとき、その逆行列も何かしらの分布に従う。その逆行列が従う分布が逆Wishart分布
- MCMCpack::riwish()関数は逆Wishart分布に従うランダム行列発生関数だが、Wishart分布行列を作成してその逆行列を計算して返している。

> riwish
function (v, S) 
{
    return(solve(rwish(v, solve(S))))
}
<environment: namespace:MCMCpack>

- $X ... N(\Sigma)$ : $X$ は分散共分散行列 $\Sigma$ が定める正規分布からの標本
- $X^TX ... Wi(\Sigma)$ : その標本のばらつき具合 $X^TX$ は $\Sigam$ が定めるWishart分布
- $Wi^{-1} (X^TX) ... Wi^{-1}(Wi(\Sigma)$ : Wishart分布の「逆関数」が逆Wishart分布なので、作用させる
- $Wi^{-1}(X^TX) ... \Sigma$ : 標本のばらつき具合 $X^TX$ を逆Wishart分布で評価すると $\Sigma$ の情報が得られる
- したがって、逆Wishart分布から得られるランダムな行列は、分散共分散行列と同じ土俵のものであるから、以下のような工夫によって、Wishart分布を「分散共分散行列を発生させる分布とみなして、その分散共分散行列に対応する多変量正規分布の広がりのばらつき」として視覚化することができる
  - 発生行列のばらつきの視覚化(ellipse()関数)
- 逆Wishart分布は分散共分散行列の共役事前分布
  - 逆Wishart分布は分散共分散行列と同じ土俵にある「分布」である
  - そこに「無情報分布」を設定したり、特定の分散共分散行列がありそうで、それを中心にある程度ばらついている、という事前予測を事前分布として入れることができる
  - 便利なのは、そのような事前分布は逆Wishart分布のパラメタとして指定できること
  - そして観測 $X$ が得られたら、 $X^TX$ とそのサンプルサイズとを情報として事後分布化するにあたり、パラメタの合成のみで行えること
Normal-Inverse-Wishart分布
- 多変量正規分布のconjugate prior
- 逆Wihart分布が分散共分散行列の共役事前分布だったが
- 多変量正規分布の期待値と分散共分散行列との両方を対象にした共役事前分布(パラメタをいじるだけで事前分布設定ができて、観測によりアップデートできる)があると便利
- それが、Normal-Inverse-Wishart分布で、正規分布と逆Wihart分布との積になっている