線形混合モデルのギブスサンプリングによるパラメータ推定① 理論編

# はじめに
機械学習の勉強方法は色々あると思いますが、私はやっぱり勉強した内容をコードで実装するまで、ちゃんと理解した感じがしないので、本を読んだ後はなるべく実装するようにしています。機械学習におけるパラメータ推定で必須のMCMCについては、メトロポリスヘイスティングは実装したことがありましたが、ギブスサンプリングは今までやったことがありませんでした。
私が所属する研究室では、遺伝子情報を用いた解析をしいていますが、そこでは線形混合モデルと呼ばれるモデルがよく出てきます。そこで今回は、この線形混合モデルに出てくるパラメータの推定をギブスサンプリングを通して行っていきたいと思います。
本投稿は、大きく二つに分けて投稿する予定です。一つ目が線形混合モデルのパラメータに対する完全事後分布を計算する理論編。もう一つが計算された完全事後分布を元にRで実装する実装編です。合わせてお読みいただければ幸いです。

[線形混合モデルのギブスサンプリングによるパラメータ推定② 実装編](https://qiita.com/KimuraSoh/items/f657847a0f354bf1c137)

### 目次
・[線形混合モデル](#線形混合モデル)
・[同時事後分布の導出](#同時事後分布の導出)
・[完全事後分布の導出](#完全事後分布の導出)

# 線形混合モデル
線形混合モデルの説明をしっかりするとそれだけで数投稿分になってしまうので、ここでは本当に簡単に線形混合モデルを復習します。まず数式ですが、次のように書き表すことができます。
```math
\displaylines{
\mathbf{y}=X\mathbf{b}+Z\mathbf{u}+\mathbf{e}\\
Var(\mathbf{u}) \sim MVN(\mathbf{0},V)\\
Var(\mathbf{e}) \sim MVN(\mathbf{0},R)
}
```
ここで、頻度主義の立場からすると、$\mathbf{b}$はある分布に従うことが仮定されないため固定効果、$\mathbf{u}$は分布に従うことが仮定されるため変量効果と呼ばれます。変量効果自身は一つの値を持たないため、頻度主義の立場から考えると、このモデルのパラメータは$\mathbf{b},V,R$となります。もちろんこれらすべてを推定することもモデル上可能ではありますが、推定するパラメータ数が多すぎてうまく推定できないことが多いため、$V$や$R$に対してはある構造を仮定してパラメータ数を減らす努力がなされます。
最もよく見かける仮定が誤差分散$R$に対する仮定で、$R=I\sigma_e^2$というように、誤差がiid(独立同一分布)であると仮定します。このようにすることでパラメータ数をひとつ($\sigma_e^2$)にすることが可能です。また、生物統計で使われるGBLUPと呼ばれるモデルでは、変量効果に対する分散共分散行列$V$に対して、$V$がゲノムから推定可能なゲノム関係行列($K$)の定数倍である($V=K\sigma_g^2$)というような仮定を置きます。このようにすることで誤差と同様、変量効果の分散共分散行列についてもパラメータ数をひとつ($\sigma_g^2$)に減らすことが可能です。つまりGBLUPで仮定されるモデルは以下の通りとなります。
```math
\displaylines{
\mathbf{y}=X\mathbf{b}+Z\mathbf{u}+\mathbf{e}\\
Var(\mathbf{u}) \sim MVN(\mathbf{0},K\sigma_g^2)\\
Var(\mathbf{e}) \sim MVN(\mathbf{0},I\sigma_e^2)
}
```
さらに簡単のため、$\lambda_g=\sigma_g^{-2},\lambda_e=\sigma_e^{-2}$と置くことで、次を得ます。
```math
\displaylines{
\mathbf{y}=X\mathbf{b}+Z\mathbf{u}+\mathbf{e}\\
Var(\mathbf{u}) \sim MVN(\mathbf{0},K/\lambda_g)\\
Var(\mathbf{e}) \sim MVN(\mathbf{0},I/\lambda_e)
}
```
今回はこのモデルの変数$\mathbf{b},\mathbf{u},\lambda_g,\lambda_e$に対して、ギブスサンプリングでパラメータ推定を行うこととします。

# 同時事後分布の導出
ギブスサンプリングを行う上でのファーストステップは同時事後分布$p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e|\mathbf{y})$を求めることです。まずは事前分布の具体形を仮定せず、事後分布が尤度と事前分布の積でどのように表せそうかを考えてみます。まずは単純なベイズの定理より次の変形が導けます。
```math
p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e|\mathbf{y})\propto  p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)
```
さらに、順々にベイズの定理を適用していくことにより次が得られます。
```math
\displaylines{
\begin{align}
p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e|\mathbf{y})&\propto p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)\\&=p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{u},\lambda_g,\lambda_e)\\&=p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{u}|\lambda_g,\lambda_e)*p(\lambda_g,\lambda_e)\\&=p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{u}|\lambda_g,\lambda_e)*p(\lambda_g|\lambda_e)*p(\lambda_e)
\end{align}
}
```
さて、まず$p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)$について考えます。本来$\mathbf{y}$の分布は、$E(\mathbf{y})=X\mathbf{b},Var(\mathbf{y})=ZGZ^T/\lambda_g+I\lambda_e$に従う正規分布となることが次の変形からわかります。
```math
\displaylines{
\begin{align}
E(\mathbf{y})&=E(X\mathbf{b}+Z\mathbf{u}+\mathbf{e})\\&=E(X\mathbf{b})+E(Z\mathbf{u})+E(\mathbf{e})\\&=X\mathbf{b}\\
Var(\mathbf{y})&=Var(X\mathbf{b}+Z\mathbf{u}+\mathbf{e})\\&=Var(X\mathbf{b})+Var(Z\mathbf{u})+Var(\mathbf{e})+2Cov(X\mathbf{b},Z\mathbf{u})+2Cov(X\mathbf{b},\mathbf{e})+2Cov(Z\mathbf{u},\mathbf{e})\\&=Var(Z\mathbf{u})+Var(\mathbf{e})\\&=ZGZ^T/\lambda_g+I\lambda_e
\end{align}
}
```
しかし、$\mathbf{u}$を固定した場合(与えられたものとした場合)、$Z\mathbf{u}$の項は固定効果と同じ扱いを受けますので、$E(\mathbf{y})=X\mathbf{b}+Z\mathbf{u},Var(\mathbf{y})=I\lambda_e$と変化します。これらの事実から、$p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)$について、次のような分布に従うことが分かります。
```math
\displaylines{
\begin{align}
p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)&=p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_e)\\&=\frac{1}{(2\pi)^{N/2}|I/\lambda_e|^{1/2}}\exp\Big(-\frac{1}{2}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(I/\lambda_e)^{-1}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})\Big)
\end{align}
}
```
ここで$N$はデータ$\mathbf{y}$の長さを表しています。

また$p(\mathbf{u}|\lambda_g,\lambda_e)$について、$\mathbf{u}$の分布は$Var(\mathbf{u}) \sim MVN(\mathbf{0},K/\lambda_g)$であることが仮定されていますので、$\mathbf{u}$の長さを$n$とすると次のように表されます。
```math
\displaylines{
\begin{align}
p(\mathbf{u}|\lambda_g,\lambda_e)&=p(\mathbf{u}|\lambda_g)\\&=\frac{1}{(2\pi)^{n/2}|K/\lambda_g|^{1/2}}\exp\Big(-\frac{1}{2}\mathbf{u}^T(K/\lambda_g)^{-1}\mathbf{u}\Big)
\end{align}
}
```
よって、同時事後分布について、
```math
\displaylines{
\begin{align}
p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e|\mathbf{y})& \propto p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e)*p(\mathbf{u}|\lambda_g,\lambda_e)*p(\lambda_g|\lambda_e)*p(\lambda_e)\\&=p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_e)*p(\mathbf{u}|\lambda_g)*p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e)*p(\lambda_g|\lambda_e)*p(\lambda_e)\\
\end{align}
}
```
where
```math
\displaylines{
\begin{align}
p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_e)*p(\mathbf{u}|\lambda_g)=\frac{1}{(2\pi)^{N/2}|I/\lambda_e|^{1/2}}\exp\Big(-\frac{1}{2}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(I/\lambda_e)^{-1}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})\Big)*\frac{1}{(2\pi)^{n/2}|K/\lambda_g|^{1/2}}\exp\Big(-\frac{1}{2}\mathbf{u}^T(K/\lambda_g)^{-1}\mathbf{u}\Big)
\end{align}
}
```
となることが分かります。次に$\mathbf{b},\lambda_g,\lambda_e$の分布を考えるのですが、基本的に、これらは事後分布と共役になるように選びます。まず$\lambda_g,\lambda_e$が考えやすいので、それらを先に考えます。
先ほど出てきた$
p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_e)*p(\mathbf{u}|\lambda_g)$において、$\lambda_g,\lambda_e$が関わっている部分のみを抜き出すと、それぞれ以下のようになります。
```math
\displaylines{
\begin{align}
&\lambda_g^{N/2}\exp\Big(-\frac{1}{2}\lambda_g(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})\Big)\\
&\lambda_e^{n/2}\exp\Big(-\frac{1}{2}\lambda_e(\mathbf{u}^TK^{-1}\mathbf{u})\Big)
\end{align}
}
```
これはガンマ分布と同じ形をしていますので、$\lambda_g,\lambda_e$について、その分布はガンマ分布を仮定すればよいことが分かります。すなわち
```math
\displaylines{
\begin{align}
p(\lambda_g|\lambda_e)=p(\lambda_g)=Gam(\lambda_g|a_g,b_g)\\
p(\lambda_e)=p(\lambda_e)=Gam(\lambda_e|a_e,b_e)
\end{align}
}
```
となります。一方の$\mathbf{b}$ですが、基本的に固定効果には特殊な事前分布をもうけないので定数を仮定することにします。すなわち
```math
\displaylines{
\begin{align}
p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e)=p(\mathbf{b})=Const
\end{align}
}
```
これで、事後分布を構成するすべての要素が出そろいました。これらをまとめると、事後分布$p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e|\mathbf{y})$は次のように書くことができそうです。

```math
\displaylines{
\begin{align}
p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e|\mathbf{y})& \propto p(\mathbf{y}|\mathbf{b},\mathbf{u},\lambda_e)*p(\mathbf{u}|\lambda_g)*p(\mathbf{b})*p(\lambda_g)*p(\lambda_e)\\&=\frac{1}{(2\pi)^{N/2}|I/\lambda_e|^{1/2}}\exp\Big(-\frac{1}{2}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(I/\lambda_e)^{-1}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})\Big)*\frac{1}{(2\pi)^{n/2}|K/\lambda_g|^{1/2}}\exp\Big(-\frac{1}{2}\mathbf{u}^T(K/\lambda_g)^{-1}\mathbf{u}\Big)*Gam(\lambda_g|a_g,b_g)*Gam(\lambda_e|a_e,b_e)
\end{align}
}
```
なお、$\mathbf{b}$に対する定数は議論に影響を及ぼさないため、省略してしまいました。ここからは、いま導出された同時分布から完全事後分布を導出していきます。

# 完全事後分布の導出
**$\mathbf{b}$に対する完全事後分布**
まずは、$\mathbf{b}$に対する完全事後分布を導出します。これはそんな難しくなく、$p(\mathbf{b},\mathbf{u},\lambda_g,\lambda_e|\mathbf{y})$から$\mathbf{b}$に関連する部分を抜き出して、正規化すれば導出できます。
$\mathbf{b}$に関連する部分を抜き出すと、
```math
\displaylines{
\begin{align}
&\exp\Big(-\frac{1}{2}\lambda_e(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})\Big)\\
= \quad &\exp\Big(-\frac{1}{2}\lambda_e(\mathbf{b}^TX^TX\mathbf{b}-2\mathbf{b}^TX^T(\mathbf{y}-Z\mathbf{u})\Big)+Const\\
= \quad &\exp\Big(-\frac{1}{2}\lambda_e\big(\mathbf{b}-(X^TX)^{-1}X^T(\mathbf{y}-Z\mathbf{u})\big)^TX^TX(\mathbf{b}-(X^TX)^{-1}X^T\big(\mathbf{y}-Z\mathbf{u})\big)\Big)+Const
\end{align}
}
```
よって$\mathbf{b}$の完全事後分布については、$p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e,\mathbf{y})= MVN(\lambda_e(X^TX)^{-1}X^T(\mathbf{y}-Z\mathbf{u}),(X^TX\lambda_e)^{-1})$とわかります。

**$\mathbf{u}$に対する完全事後分布**
次に、$\mathbf{u}$に対する完全事後分布を考えます。$\mathbf{u}$に関連する部分を抜き出すと、
```math
\displaylines{
\begin{align}
&\exp\Big(-\frac{1}{2}\lambda_e(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})\Big)*\exp\Big(-\frac{1}{2}\lambda_g(\mathbf{u}^TK^{-1}\mathbf{u})\Big)\\
= \quad &\exp\Big(-\frac{1}{2}\lambda_e(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})-\frac{1}{2}\lambda_g(\mathbf{u}^TK^{-1}\mathbf{u})\Big)\\
= \quad &\exp\Big(-\frac{1}{2}(\lambda_e\mathbf{u}^TZ^TZ\mathbf{u}-2\lambda_e\mathbf{u}^TZ^T(\mathbf{y}-X\mathbf{b})+\lambda_g\mathbf{u}^T(K)^{-1}\mathbf{u})\Big)+Const\\
= \quad &\exp\Big(-\frac{1}{2}\big(\mathbf{u}^T(Z^T\lambda_eZ+(K/\lambda_g)^{-1})\mathbf{u}-2\lambda_e\mathbf{u}^TZ^T(\mathbf{y}-X\mathbf{b})\big)\Big)+Const\\
= \quad &\exp\Big(-\frac{1}{2}\big((\mathbf{u}-(Z^T\lambda_eZ+(K/\lambda_g)^{-1})^{-1}\lambda_eZ^T(\mathbf{y}-X\mathbf{b}))^T(Z^T\lambda_eZ+(K/\lambda_g)^{-1})(\mathbf{u}-(Z^T\lambda_eZ+(K/\lambda_g)^{-1})^{-1}\lambda_eZ^T(\mathbf{y}-X\mathbf{b}))\big)\Big)+Const
\end{align}
}
```
よって$\mathbf{u}$の完全事後分布については、$p(\mathbf{u}|\mathbf{b},\lambda_g,\lambda_e,\mathbf{y})= MVN((Z^T\lambda_eZ+(K/\lambda_g)^{-1})^{-1}\lambda_eZ^T(\mathbf{y}-X\mathbf{b}),(Z^T\lambda_eZ+(K/\lambda_g)^{-1})^{-1})$とわかります。
**$\lambda_g$に対する完全事後分布**
次に、$\lambda_g$に対する完全事後分布を考えます。$\lambda_g$に関連する部分を抜き出すと、
```math
\displaylines{
\begin{align}
&\lambda_g^{N/2}\exp\Big(-\frac{1}{2}\lambda_g(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})\Big)\lambda_g^{a_g-1}\exp(-b_g\lambda_g)\\
= \quad &\lambda_g^{N/2+a_g-1}\exp\Big(-\lambda_g\big(\frac{1}{2}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})+b_g\big)\Big)
\end{align}
}
```
よって$\lambda_g$に対する完全事後分布については、$p(\lambda_g|\mathbf{b},\mathbf{u},\lambda_e,\mathbf{y})=Gam(\lambda_g|(N/2+a_g),\frac{1}{2}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})+b_g))$とわかります。

**$\lambda_e$に対する完全事後分布**
次に、$\lambda_e$に対する完全事後分布を考えます。$\lambda_e$に関連する部分を抜き出すと、
```math
\displaylines{
\begin{align}
&\lambda_e^{n/2}\exp\Big(-\frac{1}{2}\lambda_e(\mathbf{u}^TK^{-1}\mathbf{u})\Big)\lambda^{a_e-1}\exp(-b_e\lambda_e)\\
= \quad &\lambda_e^{n/2+a_e-1}\exp\Big(-\lambda_e\big(\frac{1}{2}(\mathbf{u}^TK^{-1}\mathbf{u})+b_e\big)\Big)
\end{align}
}
```
よって$\lambda_e$に対する完全事後分布については、$p(\lambda_e|\mathbf{b},\mathbf{u},\lambda_g,\mathbf{y})=Gam(\lambda_e|(n/2+a_e),\frac{1}{2}(\mathbf{u}^TK^{-1}\mathbf{u})+b_e)$とわかります。


# まとめ
これまでの計算を通して、各変数の完全事後分布が以下のようにもとまりました。

```math
\displaylines{
\begin{align}
&p(\mathbf{b}|\mathbf{u},\lambda_g,\lambda_e,\mathbf{y})= MVN(\lambda_e(X^TX)^{-1}X^T(\mathbf{y}-Z\mathbf{u}),(X^TX\lambda_e)^{-1})\\
&p(\mathbf{u}|\mathbf{b},\lambda_g,\lambda_e,\mathbf{y})= MVN((Z^T\lambda_eZ+(K/\lambda_g)^{-1})^{-1}\lambda_eZ^T(\mathbf{y}-X\mathbf{b}),(Z^T\lambda_eZ+(K/\lambda_g)^{-1})^{-1})\\
&p(\lambda_g|\mathbf{b},\mathbf{u},\lambda_e,\mathbf{y})=Gam(\lambda_g|(N/2+a_g),\frac{1}{2}(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})^T(\mathbf{y}-X\mathbf{b}-Z\mathbf{u})+b_g))\\
&p(\lambda_e|\mathbf{b},\mathbf{u},\lambda_g,\mathbf{y})=Gam(\lambda_e|(n/2+a_e),\frac{1}{2}(\mathbf{u}^TK^{-1}\mathbf{u})+b_e)
\end{align}
}
```
ギブスサンプリングでは、これらの完全事後分布を使い、次のようにサンプリングを行います。

```math
\displaylines{
\begin{align}
&\mathbf{b}_i\sim p(\mathbf{b}|\mathbf{u}_{i-1},{\lambda_g}_{i-1},{\lambda_e}_{i-1},\mathbf{y})\\
&\mathbf{u}_i\sim p(\mathbf{u}|\mathbf{b}_i,{\lambda_g}_{i-1},{\lambda_e}_{i-1},\mathbf{y})\\
&{\lambda_g}_i\sim p(\lambda_g|\mathbf{b}_i,\mathbf{u}_i,{\lambda_e}_{i-1},\mathbf{y})\\
&{\lambda_e}_i\sim p(\lambda_e|\mathbf{b}_i,\mathbf{u}_i,{\lambda_g}_i,\mathbf{y})\\

\end{align}
}
```
次回の投稿では実際にこれらを実装し、うまく動くことを確認していきます。なお、ここまでの計算はかなり適当に進んできたため、もしかすると計算過程に間違いがあるかもしれません。お気づきの際はご連絡いただけますと幸いです。

**次の投稿**
[線形混合モデルのギブスサンプリングによるパラメータ推定② 実装編](https://qiita.com/KimuraSoh/items/f657847a0f354bf1c137)