BLUPの導出 ~平均二乗誤差から~

## はじめに
本投稿は、佐々木義之さんの「変量効果の推定とBLUP法」の第1章(p20-25)で行われているBLUPの導出に関して、誤差分散の最小化によりではなく、平均二乗誤差の最小化により、導出を試みたものになります。BLUPの定義自体は誤差分散の最小化ですが、本教科書では、二乗和による最小二乗推定量の導出、平均二乗誤差の最小化によるBLPの導出、という流れを分でBLUPを導出しておりますので、平均二乗誤差の最小化を通したBLUPの導出を行った方がすっきりしています。BLUPでは不偏性を仮定しておりますので、誤差分散=平均二乗誤差であり、平均二乗誤差からBLUPが求められるのは当然なのですが、少し計算がややこしいので、本投稿でちゃんと導出を行います。

## 誤差分散の最小化によるBLUPの導出
まず、誤差分散の最小化によるBLUPの導出を復習しておきます。
問題設定としては以下の通りです。

①$\mathbf{y}$は$X\mathbf{b}$の不偏推定量の一つ ($E[\mathbf{y}]=X\mathbf{b}$)
②$\mathbf{y}$の分散は$V$ ($Var(\mathbf{y})=V$)
この条件を満たす$\mathbf{y}$とベクトル$\mathbf{a}$の内積$\mathbf{a}^T\mathbf{y}$を用いて推定可能関数$\mathbf{k}^T\mathbf{b}$のBLUPを求めたい。

BLUPは①線形で②不偏性を持ち③誤差分散を最小化する推定量であり、①は$\mathbf{a}^T\mathbf{y}$を仮定した時点で成り立っています。また、②については、$\mathbf{k}^T=\mathbf{a}^TX$という条件と等価です。したがって、$Var(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})$を$\mathbf{k}^T=\mathbf{a}^TX$という条件を満たすように、最小化する$\mathbf{a}$が、BLUPであることが直ちにわかります。これにはラグランジュの未定乗数法を用いればよく、最小化すべき関数$f(\mathbf{a},\mathbf{t})$は次のように求まります。

```math
\begin{align}
f(\mathbf{a},\mathbf{t})&=Var(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})\\
&=\mathbf{a}^TVar(\mathbf{y})\mathbf{a}+Var(\mathbf{k}^T\mathbf{b})-2Cov(\mathbf{a}^T\mathbf{y},\mathbf{k}^T\mathbf{b})+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})\\
&=\mathbf{a}^TV\mathbf{a}-2\mathbf{a}^TCov(\mathbf{y},\mathbf{k}^T\mathbf{b})+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})
\end{align}
```
最後の式への式変形では、$\mathbf{a},\mathbf{t}$どちらとも関係ない$Var(\mathbf{k}^T\mathbf{b})$の項を消去しました。

ここからの計算は、[別の投稿](https://qiita.com/KimuraSoh/items/f9a7212a22ddf5fcb655)に詳しく記述しているので、ここでは省略しますが、最後の式の形にまで持っていければ、BLUPが求められることを覚えておいてください。

## 平均二乗誤差の最小化によるBLUPの導出
さて、先ほど用いた誤差分散$Var(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})$ですが、分散の一般的な関係式より次のように変形することができます。
```math
\begin{align}
Var(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})
&=E[(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})^2]-(E[\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b}])^2
\end{align}
```
ここで、右辺第二項目の$(E[\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b}])^2$に着目します。期待値についての公式を使い変形すると以下を得ます。
```math
\begin{align}
(E[\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b}])^2
&=(E[\mathbf{a}^T\mathbf{y}]-E[\mathbf{k}^T\mathbf{b}])^2\\
&=(\mathbf{k}^T\mathbf{b}-E[\mathbf{k}^T\mathbf{b}])^2\\
\end{align}
```
これは$\mathbf{a},\mathbf{t}$どちらとも関係ありませんので、その意味で定数と見なせます。したがって、これを$C$とすると、
```math
\begin{align}
Var(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})
&=E[(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})^2]+C
\end{align}
```
これは誤差分散と平均二乗誤差が($\mathbf{a},\mathbf{t}$に対する微分の意味で)、同等であることを意味しています。そこで、平均二乗誤差に着目し、その最小化によりBLUPが求まるかを考えてみます。

考えるべき関数$f(\mathbf{a},\mathbf{t})$は次の通りです。
```math
\begin{align}
f(\mathbf{a},\mathbf{t})=E[(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})^2]+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})
\end{align}
```
これを変形して、
```math
\begin{align}
f(\mathbf{a},\mathbf{t})&=\mathbf{a}^TV\mathbf{a}-2\mathbf{a}^TCov(\mathbf{y},\mathbf{k}^T\mathbf{b})+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})
\end{align}
```
とすることを目指します。しかし、それは比較的簡単で、以下のようにすぐ求まります。
```math
\begin{align}
f(\mathbf{a},\mathbf{t})&=E[(\mathbf{a}^T\mathbf{y}-\mathbf{k}^T\mathbf{b})^2]+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})\\
&=E[\mathbf{a}^T\mathbf{y}\mathbf{y}^T\mathbf{a}+\mathbf{k}^T\mathbf{b}\mathbf{b}^T\mathbf{k}+2\mathbf{a}^T\mathbf{y}\mathbf{b}^T\mathbf{k}]+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})\\
&=\mathbf{a}^TV\mathbf{a}+\mathbf{a}^T\mathbf{b}^TX^TX\mathbf{b}\mathbf{a}+E[\mathbf{k}^T\mathbf{b}\mathbf{b}^T\mathbf{k}]+2Cov[\mathbf{a}^T\mathbf{y},\mathbf{b}^T\mathbf{k}]+2E[\mathbf{a}^T\mathbf{y}]E[\mathbf{b}^T\mathbf{k}]+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})\\
&=\mathbf{a}^TV\mathbf{a}+\mathbf{a}^T\mathbf{b}^TX^TX\mathbf{b}\mathbf{a}+E[\mathbf{k}^T\mathbf{b}\mathbf{b}^T\mathbf{k}]+2\mathbf{a}^TCov[\mathbf{y},\mathbf{b}^T\mathbf{k}]+2\mathbf{b}^T\mathbf{k}E[\mathbf{b}^T\mathbf{k}]+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})\\
&=\mathbf{a}^TV\mathbf{a}+2\mathbf{a}^TCov[\mathbf{y},\mathbf{b}^T\mathbf{k}]+2\mathbf{t}^t(\mathbf{k}-X^T\mathbf{a})\\
\end{align}
```
以上より、平均二乗誤差の最小化からもBLUPが求められることが分かりました。