# 黒木玄(@genkuroki)氏によるガウス過程とベイズ法に関する説明

https://twitter.com/genkuroki/status/1476123484069965824 からの一連のスレッドによる。しかしところどころ@gm3d2による勝手な文言を付け加えていたり言い回しを変えているので注意

## 普遍的に用いられるアプローチ

1. モデルのパラメータの値をなんらかの「尤度」最大化によって調節する。
2. モデルの中でデータと同じ値が生じたという条件で定義される条件付き確率分布を考える。

いわゆるベイズ統計、ガウシアン回帰などの各種法は多くの場合これの特殊な事例として見ることができる。
(前提として、あくまでモデルに依存した議論であることには注意する)

## 出発点のモデル
ここでは変数$x_i$のそれぞれに対して値$y_i$が観測される確率の密度が以下で与えられるとする
<!-- $$ Z(y_1,…,y_n|x_1,…,x_n, heta), ~n=0,1,2,… $$ -->

$$ Z(y_1,\ldots,y_n | x_1,\ldots,x_n, \theta) $$

データ $(x_1, y_1),\ldots,(x_n, y_n)$ が得られたとき、1. の方法に沿ってモデルの尤度関数
$$ \theta \mapsto Z(y_1,\ldots,y_n | x_1,\ldots,x_n, \theta) \tag{1}\label{eq1}$$
を最大化するようにパラメータ$\theta$の値を決める。最尤法もハイパーパラメータの調節もこれの特別な場合に過ぎない

2. の方針に従うと、既知のデータ$(x_1, y_1),\ldots,(x_n, y_n)$とその次に与えられた$x_\ast$から得られるデータ$y_\ast$の予測分布は、

**モデル内で既知のデータの値と同じ値が生成されたという条件のもとでの、次に得られるデータの値に関する条件付き確率分布**

として定義され、次の形になる

$$Z(y_\ast|x_\ast,(x_1, y_1),\ldots,(x_n,y_n), \theta)
= Z(y_1,\ldots,y_n,y_\ast|x_1,\ldots,x_n,x_\ast,\theta) / Z(y_1,\ldots,y_n|x_1,\dots,x_n, \theta)$$

**これ自体は条件付確率の定義と各$(x_i, y_i)$の独立性から簡単に得られる**

## ベイズ統計における事前分布や事後分布はどうなったのか？
上記の特別なケースとして、モデルがそれらを使って記述される場合に出てくる

### 元のモデルが事前分布に基づいて与えられている例

$$ Z(y_1,\ldots,y_n | x_1,\ldots,x_n,\theta)
=\int p(y_1|x_1,\theta,\eta)\ldots p(y_n|x_n,\theta,\eta) \phi(\eta|\theta)d\eta $$

(i.i.d. の場合)

$p(y_i|x_i,\theta,\eta)$: パラメータ付き確率密度関数、$\phi(\eta|\theta)$: 事前分布

モデル内で実際に登場しているのは$y_1,…,y_n,\eta$の同時確率密度関数$p(y_1|x_1,\theta,\eta)\ldots p(y_n|x_n,\theta,\eta) \phi(\eta|\theta)$であることに注意

$x_1,\ldots,x_n$および$\theta$は外部から見えるパラメータ扱い

この状況で既知のデータと同じ値がモデル内部で生成されたという条件のもとでの内部パラメータηの条件付き確率分布の密度関数は

$$ \phi(η|(x_1,y_1),\ldots,(x_1,y_n),\theta)
= p(y_1|x_1,\theta,\eta)…p(y_n|x_n,\theta,η)φ(η|\theta) / Z(y_1\dots,y_n|x_1,\ldots,x_n,\theta)$$

となり、これがいわゆる事後分布。これを使うと$y_*$の事後予測分布は

$$ Z(y_\ast|x_\ast,(x_1,y_1)\ldots,(x_n,y_n), \theta)
=\int p(y_*|\theta,η) \phi(η|(x_1,y_1),\ldots,(x_1,y_n),\theta) d\eta$$

とも書け、いわゆるベイズ法での事後予測分布の形が得られる。**つまり個別にベイズ法における事後分布を定義として導入する必要はない**

## この視点での最尤法

元のモデルが事前分布を用いない形の
$$ Z(y_1,\ldots,y_n|x_1,\ldots,x_n,\theta)=p(y_1|x_1,\theta)…p(y_n|x_n,\theta) $$
で与えられており、パラメータ$\theta$を「尤度」$Z(y_1,\ldots,y_n|x_1,\ldots,x_n,\theta)$の最大化によって決めようとする場合に相当

このとき事後予測分布$Z(y_\ast|x_\ast,(x_1,y_1),\ldots,(x_n,y_n),\theta)$を上の一般論に従って計算すると簡単な形になり、
$$ Z(y_1,\ldots,y_n,y_\ast|x_1,\ldots,x_n,x_\ast,\theta)/Z(y_1,\ldots,y_n|x_1,\ldots,x_n,\theta) = p(y_\ast|x_\ast,\theta) $$
と一変数についてのモデルによる確率密度関数そのものになる。これを$\theta$の関数と見て最大にする$\hat{\theta}$を採用するとよく知られた最尤法での予測分布になる

以上のように**モデル内条件付き確率分布としての予測分布**に注目すればベイズ法も最尤法も共通の枠組みの特殊な場合として見ることができる

## 補足1
「尤度最大化」には「正則化された(罰則付きの)尤度最大化」という一般化がある。$\eqref{eq1}$を適当な正値関数（罰則関数）$\psi(\theta)$を用いて
$$\theta \mapsto Z(y_1,\ldots,y_n | x_1,\ldots,x_n,\theta) \psi(\theta)$$と拡張することに相当

特に正値関数$\psi(\theta)$が$\theta$に関する確率密度関数になっていれば事前分布とみなすことができ、
$$ Z(y_1,\ldots,y_n|x_1,\ldots,x_n)=\int Z(y_1,\dots,y_n|x_1,\ldots,x_n,\theta)\psi(\theta)d\theta $$
で外部に調節するパラメータが見えていないモデルを作れ、いわゆる「フルベイズ」と言われるものになる

ただし、正値関数と確率密度関数では、座標変換での変換性が違うので、その点に注意してしないと、正則化法(事後確率最大化法)とベイズ法の関係を座標変換時に誤解する危険性がある

ベイズ法に使う事前分布を尤度関数の正則化因子扱いする場合には、座標変換時に確率密度関数扱いするのではなく、単なる正値関数扱いする必要がある

(渡辺澄夫『ベイズ統計の理論と方法』も参照)

パラメータをどのように調節するかについては方針が無限にあるのでそれに応じて$\psi(\theta)$の選択肢も無限にある

パラメータの個数がデータの個数に対して小さくない場合、単純な最尤法だとパラメータの値が大きく暴れやすい（overfitting)。これを$\psi(\theta)$による$\theta$に対する重み付けによって制御するのが多くの場合に行われる正則化

実際には正則化因子$\psi(\theta)$もパラメータ$\lambda$(ハイパーパラメータ)を持つ$\psi(\theta,\lambda)$の形をしていて、データを用いてそのパラメータ$\lambda$を調節する方法(クロスバリデーションなど)を決めておくのが普通

### 古典的で有名な正則化法の例
#### Stein推定

nbviewer.org/gist/genkuroki…

を参照。

## 補足2
分散は同じだが平均値は異なるかもしれない独立なn個の正規分布で生成された乱数から、n個の正規分布の未知の平均値を推定するときには、最尤法は平均二乗誤差的に損な方法になる

## 補足3
### Gauss過程回帰

$$ Z(y_1,\ldots,y_n|x_1,\ldots,x_n,\theta,\sigma^2)
=(2π \det(\Sigma))⁻¹ᐟ²
×\exp(-(1/2)(y-\mu(x))^T \Sigma⁻¹ (y-\mu(x))) $$

ここで、

$y = \begin{pmatrix}y_1\\ \vdots\\y_n\end{pmatrix}$
$\mu(x) = \begin{pmatrix}\mu(x_1)\\ \vdots \\ \mu(x_n)\end{pmatrix}$
$\Sigma = (\Sigma_{ij}) = (k(x_i, x_j|\theta)) + \sigma^2 I)$ (正方行列)

は関数$\mu$と$k$を与えると決まる。これに一般的な方針1.、2.を適用するとGauss過程回帰になる。

- 多くのモデルの構成はこのスレッドの記号での$Z$の記述に集約される
- 例えば、Stanで計算するようなモデルはこの範疇に全部入っている
- $Z$を正確に記述して、その直観的意味についてよく考えることは有益
- $Z$は統計力学での分配関数に似ている

コンピュータによる具体的計算をしたければ$\mu(x)=0$、$k(x,x^\prime|\theta) = a * exp(-(x - x^\prime)²/(2b²))$、$\theta = (a, b)$として試してみるとよい

- やっていることは単に「有限個の与えられた点$x_i$たちにおける値$y_i$たちが多変量正規分布に従う」と考え、一般的な方針1.、2.を適用しているだけ
- $Z$が多変量正規分布の密度関数そのものであることに注意し、上記の具体的なGaussカーネル$k$の場合に対応する多変量正規分布の乱数のサンプルをプロットしてみると雰囲気がよくわかる