はじめに

このノートは「StanとRでベイズ統計モデリング」の内容を写経することで、ベイズ統計への理解を深めていくために作成している。

基本的には気になった部分を写経しながら、ところどころ自分用の補足をメモすることで、「StanとRでベイズ統計モデリング」を読み進めるための自分用の補足資料になることを目指す。私の解釈がおかしく、メモが誤っている場合があるので注意。

今回は第10章「収束しない場合の対処法」の後半、10.2.4の分散共分散行列のチャプターから写経していく。

10.2.4 分散共分散行列

チャプター8の階層モデルを多変量正規分布、共分散行列を使って拡張する。このモデルのメカニズムとしては、$a[1],a[2],a[3],a[4]$を「すべての会社で共通の全体平均」と「会社差を表す項」に分解する。会社差は平均0標準偏差$\sigma_{a}$の正規分布から生成されると考える。$b[1],b[2],b[3],b[4]$も同様で、年齢に伴う昇給額は「すべての会社で共通の全体平均」と「会社差を表す項」に分解する。

このような制約を入れることで、新卒時点の年収の会社差のばらつきは$\sigma_{a}$ぐらいと解釈できるようになる。$\sigma_{a}$には無情報事前分布を仮定する。このように階層的に事前分布が設定されているため階層モデルと呼ばれる。

モデル8-4

\[ \begin{eqnarray} Y[n] &\sim& Normal(a[KID[n]] + b[KID[n]]X[n], \sigma_{Y}) \\ a[k] &\sim& Normal(a_{全体平均}, \sigma_{a}) \\ b[k] &\sim& Normal(b_{全体平均}, \sigma_{b}) \\ \end{eqnarray} \]

上記モデルでは、$a[k],b[k]$は独立した正規分布に従うと仮定されている。ただ、モデルのメカニズムとして、新卒年収が高いと年齢による昇給額も高くなる高くなりやすい場合もある。$a[k]$が高いと$b[k]$も高い、$a[k]$が高いと$b[k]$が低い、となるようような相関関係をモデル化したい。このようなケースで多変量正規分布が役立つ。

モデル10-5

一般に、多変量正規分布は下記の通り表現される。この分布は$\overrightarrow{y}$の各要素$y_{1},...,y_{n}$が従う同時分布になっている。$\Sigma$は共分散行列で、$|\Sigma|$は行列式である。共分散行列には無情報事前分布を利用する。

\[ MultiNormal(\overrightarrow{y}|\overrightarrow{\mu},\Sigma) = \frac{1}{(2\pi)^{\frac{K}{2}}}\frac{1}{\sqrt{|\Sigma|}} \exp \left( -\frac{1}{2} (\overrightarrow{y}-\overrightarrow{\mu})\Sigma^{-1}(\overrightarrow{y}-\overrightarrow{\mu}) \right) \]

$\Sigma$の対角成分以外が0である場合、互いに独立な$K$個の正規分布から生成されたのと同じである。2つの要素しかない場合、$\overrightarrow{y}=(a,b)^{T}$を考えるとわかり良い。$\mu_{a},\sigma_{a}$を$a$の平均、標準偏差、$\mu_{b},\sigma_{b}$を$b$の平均、標準偏差とする。$a,b$の相関を$\rho$とする。

\[ \overrightarrow{\mu} = \left( \begin{array}{c} \mu_{a} \\ \mu_{b} \\ \end{array} \right) , \Sigma = \begin{pmatrix} \sigma_{a}^{2} & \sigma_{a}\sigma_{b}\rho \\ \sigma_{a}\sigma_{b}\rho & \sigma_{b}^{2} \end{pmatrix} \]

Stanのモデル式は下記の通り。

data {
  int N;
  int K;
  vector[N] X;
  vector[N] Y;
  int<lower=1, upper=K> KID[N];
}

parameters {
  vector[2] ab[K];
  vector[2] ab0;
  cov_matrix[2] cov;
  real<lower=0> s_Y;
}

transformed parameters {
  vector[K] a;
  vector[K] b;
  for (k in 1:K) {
    a[k] = ab[k,1];
    b[k] = ab[k,2];
  }
}

model {
  ab ~ multi_normal(ab0, cov);
  Y ~ normal(a[KID] + b[KID] .* X, s_Y);
}
// forloop版
// model {
//   for (k in 1:K) {
//     ab[k] ~ multi_normal(ab0, cov);
//   }
//   for (n in 1:N) {
//     Y[n] ~ normal(a[KID[n]] + b[KID[n]] * X[n], s_Y);
//   }
// }

ただ、このモデルは収束しない。モデルが複雑でパラメタの制限が少ないのに対し、データが十分ではないことが原因。このような場合、分散共分散行列に弱情報事前分布を設定する。このようなケースでは逆ウィシャート分布が使われることが多い模様。逆ウィシャート分布は逆ガンマ分布の多変量版で、共分散行列の共役事前分布でもある。逆Wishart分布は、Wishart分布に従う行列の逆行列が従う分布です。行列のサイズを$K$としたとき、対称正定値行列$\Sigma$の逆ウィシャート分布は、自由度パラメタ$\nu \in (K-1, \infty)$と対称正定値行列のパラメタ$S \in \mathbb{ R }^{K×K}$を使って下記のように表現される。

\[ p(\mathbf {\Sigma} | \nu, S) = \frac{1}{2^\frac{\nu K}{2}} \ \frac{1}{\Gamma_K \! \left( \frac{\nu}{2} \right)} \ \left| S \right|^\frac{\nu}{2} \ \left| \mathbf {\Sigma} \right|^{-\frac{\nu + K + 1}{2}} \ \exp \! \left( - \frac{1}{2} \ \text{tr} (S \mathbf {\Sigma}^{-1}) \right) \]

ただ、逆ウィシャート分布はパラメタの設定が難しく、解決策として、標準偏差と相関から分散共分散行列を構成する方法がある。こうすることで、$\sigma_{a},\sigma_{b}$に事前分布を設定しやすく、$\sigma_{a},\sigma_{b},\rho$から分散共分散行列$\Sigma$を構成する。

モデル10-6

\[ \begin{eqnarray} Y[n] &\sim& Normal(a[KID[n]] + b[KID[n]]X[n], \sigma_{Y}) \\ \begin{pmatrix} a[k] \\ b[k] \\ \end{pmatrix} &\sim& MultiNormal \left( \begin{array}{c} \left( \begin{array}{c} a_{全体平均} \\ b_{全体平均} \\ \end{array} \right) \end{array}, \begin{pmatrix} \sigma_{a}^{2} & \sigma_{a}\sigma_{b}\rho \\ \sigma_{a}\sigma_{b}\rho & \sigma_{b}^{2} \end{pmatrix} \right) \\ a_{全体平均} &\sim& Normal(400,200) \\ b_{全体平均} &\sim& Normal(15,15) \\ \sigma_{a} &\sim& Student\_t^{+}(4,0,200) \\ \sigma_{b} &\sim& Student\_t^{+}(4,0,20) \\ \end{eqnarray} \]

事前分布を利用するにあたり、下記の通り仮定を置いている。

新卒の基本年収の全体平均$a_{全体平均}$: 200-600万円の範囲だろう
年齢に伴う昇給額の全体平均$b_{全体平均}$: 0-30万円の範囲だろう
新卒の基本年収の会社差のばらつき$\sigma_{a}$: 高々200万円くらいだろう
年齢に伴う昇給額の会社差のばらつき$\sigma_{b}$: 高々20万円くらいだろう
$\sigma_{a},\sigma_{b}$の関係を表す$\rho$には無情報事前分布を利用

data {
  int N;
  int K;
  vector[N] X;
  vector[N] Y;
  int<lower=1, upper=K> KID[N];
}

parameters {
  vector[2] ab[K];
  vector[2] ab0;
  real<lower=0> s_a;
  real<lower=0> s_b;
  real<lower=-1, upper=1> rho;
  real<lower=0> s_Y;
}

transformed parameters {
  vector[K] a;
  vector[K] b;
  matrix[2,2] cov;
  for (k in 1:K) {
    a[k] = ab[k,1];
    b[k] = ab[k,2];
  }
  // 共分散分散行列
  cov[1,1] = square(s_a); cov[1,2] = s_a*s_b*rho;
  cov[2,1] = s_a*s_b*rho; cov[2,2] = square(s_b);
}

model {
  ab0[1] ~ normal(400, 200);
  ab0[2] ~ normal(15, 15);
  s_a ~ student_t(4, 0, 200);
  s_b ~ student_t(4, 0, 20);
  ab ~ multi_normal(ab0, cov);
  Y ~ normal(a[KID] + b[KID] .* X, s_Y);
}

library(dplyr)
library(rstan)
library(ggplot2)

options(max.print = 999999)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

library(rstan)

d <- read.csv('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap08/input/data-salary-2.txt')
N <- nrow(d)
K <- 4
data <- list(N=N, K=K, X=d$X, Y=d$Y, KID=d$KID)

モデルは下記の通りである。

data {
  int N;
  int K;
  vector[N] X;
  vector[N] Y;
  int<lower=1, upper=K> KID[N];
}

parameters {
  vector[2] ab[K];
  vector[2] ab0;
  real<lower=0> s_a;
  real<lower=0> s_b;
  real<lower=-1, upper=1> rho;
  real<lower=0> s_Y;
}

transformed parameters {
  vector[K] a;
  vector[K] b;
  matrix[2,2] cov;
  for (k in 1:K) {
    a[k] = ab[k,1];
    b[k] = ab[k,2];
  }
  // Σを構成
  cov[1,1] = square(s_a); cov[1,2] = s_a*s_b*rho;
  cov[2,1] = s_a*s_b*rho; cov[2,2] = square(s_b);
}

model {
  // 事前分布
  ab0[1] ~ normal(400, 200);
  ab0[2] ~ normal(15, 15);
  s_a ~ student_t(4, 0, 200);
  s_b ~ student_t(4, 0, 20);
  ab ~ multi_normal(ab0, cov);

  // `.*`は要素ごとの積
  Y ~ normal(a[KID] + b[KID] .* X, s_Y);
}
// forloop版
// model {
//   s_a ~ student_t(4, 0, 200);
//   s_b ~ student_t(4, 0, 20);
//   ab0[1] ~ normal(400, 200);
//   ab0[2] ~ normal(15, 15);
// 
//   for (k in 1:K) {
//     ab[k] ~ multi_normal(ab0, cov);
//   }
// 
//   for (n in 1:N) {
//     Y[n] ~ normal(a[KID[n]] + b[KID[n]] * X[n], s_Y);
//   }
// }

ここでは、stan_model()関数で最初にコンパイルしておいてから、

model106 <- stan_model('note_ahirubayes13-106.stan')

sampling()関数でサンプリングする。

fit <- sampling(object = model106, data = data, seed = 1989)

推定結果は下記の通り。

print(fit, prob = c(0.025, 0.5, 0.975), digits_summary = 1)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##             mean se_mean      sd    2.5%    50%   97.5% n_eff Rhat
## ab[1,1]    383.3     0.4    15.0   353.8  383.1   413.2  1814    1
## ab[1,2]      7.7     0.0     0.9     5.9    7.7     9.6  1736    1
## ab[2,1]    334.9     0.4    18.0   300.6  334.6   370.8  2394    1
## ab[2,2]     19.4     0.0     1.3    16.7   19.4    22.0  2424    1
## ab[3,1]    330.2     0.9    32.8   262.8  332.0   390.5  1448    1
## ab[3,2]     11.7     0.0     1.6     8.7   11.6    15.0  1577    1
## ab[4,1]    438.9     3.5   110.7   298.7  406.4   713.8  1028    1
## ab[4,2]     11.5     0.1     4.5     0.3   12.7    17.4  1037    1
## ab0[1]     370.1     1.3    57.5   264.5  363.8   507.2  2054    1
## ab0[2]      12.7     0.1     4.3     4.2   12.7    22.1  1721    1
## s_a         91.3     2.3    74.0    12.5   70.2   282.5  1042    1
## s_b          8.4     0.1     4.6     3.3    7.3    20.8  1781    1
## rho         -0.3     0.0     0.5    -0.9   -0.4     0.7  2057    1
## s_Y         28.8     0.1     3.9    22.3   28.5    37.7  1340    1
## a[1]       383.3     0.4    15.0   353.8  383.1   413.2  1814    1
## a[2]       334.9     0.4    18.0   300.6  334.6   370.8  2394    1
## a[3]       330.2     0.9    32.8   262.8  332.0   390.5  1448    1
## a[4]       438.9     3.5   110.7   298.7  406.4   713.8  1028    1
## b[1]         7.7     0.0     0.9     5.9    7.7     9.6  1736    1
## b[2]        19.4     0.0     1.3    16.7   19.4    22.0  2424    1
## b[3]        11.7     0.0     1.6     8.7   11.6    15.0  1577    1
## b[4]        11.5     0.1     4.5     0.3   12.7    17.4  1037    1
## cov[1,1] 13810.8   699.1 27249.6   157.2 4927.4 79822.1  1519    1
## cov[1,2]  -287.6    17.0   804.7 -2175.2 -129.2   807.9  2245    1
## cov[2,1]  -287.6    17.0   804.7 -2175.2 -129.2   807.9  2245    1
## cov[2,2]    92.5     3.0   130.2    10.7   53.3   431.3  1913    1
## lp__      -172.9     0.1     3.2  -180.0 -172.6  -167.4  1079    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Jan  3 16:37:53 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

この方法では、分散共分散行列のサイズが大きくなると、事前分布の設定が煩雑になる。その解決策として、LKJ相関分布を利用する方法がある。詳しくはStanマニュアルを参照。相関行列に対してLKJ相関分布で事前分布を設定し、行列演算から分散共分散行列を作る。

モデル10-7

$$ \[\begin{eqnarray} Y[n] &\sim& Normal(a[KID[n]] + b[KID[n]]X[n], \sigma_{Y}) \\ \begin{pmatrix} a[k] \\ b[k] \\ \end{pmatrix} &\sim& MultiNormal\_Cholesky \left( \begin{array}{c} \left( \begin{array}{c} a_{全体平均} \\ b_{全体平均} \\ \end{array} \right) \end{array}, \Sigma_{chol} \right) \\ \Sigma_{chol} &\sim& \begin{pmatrix} \sigma_{a} & 0 \\ 0 & \sigma_{b} \end{pmatrix} \Omega_{chol} \\ a_{全体平均} &\sim& Normal(400,200) \\ b_{全体平均} &\sim& Normal(15,15) \\ \sigma_{a} &\sim& Student\_t^{+}(4,0,200) \\ \sigma_{b} &\sim& Student\_t^{+}(4,0,20) \\ \Omega_{chol} &\sim& LKJcorr\_Cholesky(\nu) \end{eqnarray}\] $$

$\Sigma_{chol}$は分散共分散行列もコレスキー因子を表し、$\Omega_{chol}$は相関行列のコレスキー因子を表す。MultiNormal_Cholesky分布は、コレスキー因子を引数にとる多変量正規分布。LKJcorr_Cholesky分布は相関行列のコレスキー因子を生成する分布で、$\Omega_{chol}$の事前分布である。$\nu$はLKJcorr_Cholesky分布のshapeパラメタで1以上の値をとる。1だと一様分布に相当する無情報事前分布となり。2や4などであれば制約の強い相関行列が生成される。