PRML復活の呪文 part7 (2.3.7 - 2.5.2)

TL;DR

• スチューデントのt分布は見た目はガウス分布と似ているが外れ値の影響を受けにくい。
• （単一の）ガウス分布などパラメトリックな分布では表現が難しい分布は混合ガウス分布やノンパラメトリックなアプローチで表現する方法がある

2.3.7 スチューデントのt分布

平均$\mu$, 分散$\tau^{-1}$なる1変数のガウス分布を$p(x | \mu, \tau^{-1} )$とし、分散の事前分布にガンマ分布$p(\tau | a, b)$を使う。
確率のprouct ruleより$ p(x | \mu, \tau^{-1} ) p(\tau | a, b) = p(x, \tau | \mu, a, b )$なので$\tau$を消すには$\tau$について積分（周辺化）すればよい。すなわち、

$$
p(x | \mu, a, b ) = \int_0^{\infty} p(x | \mu, \tau^{-1} ) p(\tau | a, b) d\tau
$$

これを計算して得られるのがスチューデントのt分布と呼ばれる以下の式（慣例により、$\nu = 2a, \lambda = a/b$とおく）

$$
\text{St}(x | \mu, a, b) = \frac{ \Gamma(\nu/2 + 1/2) }{ \Gamma(\nu/2) }
\left( \frac{\lambda}{\pi \nu} \right)^{1/2}
\left[ 1 + \frac{ \lambda (x - \mu)^2 }{ \nu } \right]^{- \nu /2 - 1/2} \tag{2.159}
$$

イメージとしては分散0のガウス分布から分散$\infty$のガウス分布までを全部足し合わせたものの平均。

下図で緑線で表されているガウス分布より、青線や赤線で表されているスチューデントのt分布は分布の「すそ」が長い、すなわち山になっている平均から大きく離れた値をとる確率もある程度あるとするような分布となっている。

そのため、下図のように外れ値があるデータ（右側の離れた位置にあるヒストグラム）があっても影響を受けにくいという性質がある。この性質を頑健性 (robustness) とよぶ¹。

2.3.8 周期変数

省略。あんまり使うケースないし、この節の知識は今後のテキストで使わないので。

2.3.9 混合ガウス分布

ガウス分布は単峰性、つまり分布に山が１つしかない。なので、下図のように分布のかたまりが２つあるようなデータを（単一の）ガウス分布でうまくフィッティングすることが出来ない。実際、単一のガウス分布でフィッティングすると図左側のようにかたまりとかたまりの間のデータの密度が薄いところにガウス分布の山がきて、ここが最もデータが生じる確率が高いと推定している。

そこで、平均・分散が異なる$K$個のガウス分布を重ね合わせる混合ガウス分布 (mixture of Gaussians) について考える。$k$番目のガウス分布の重みに相当する係数（混合係数）を$\pi_k$として

$$
p( \mathbf{x} ) = \sum_{k=1}^K \pi_k N( \mathbf{x} | \mathbf{\mu} _ k, \Sigma_k )
\tag{2.188}
$$

重ね合わせた後の混合ガウス分布が確率、つまりxについて積分すると1になるように混合係数は

$$
\sum_{k=1}^K \pi_k = 1, 0 \leq \pi_k \leq 1 \tag{2.189, 2.190}
$$

とする。

ここから少し見方を変える。$\pi_k$を$k$番目のガウス分布が選択される事前確率、$k$番目のガウス分布$N( \mathbf{x} | \mathbf{\mu} _ k, \Sigma_k )$を$p( \mathbf{x} |k)$とおこう。
すると、ベイズの定理により

$$
p(k | \mathbf{x} ) = \frac{ p(k) p(\mathbf{x} | k) }{ \sum_l p(l) p(\mathbf{x} | l)}
\tag{2.192}
$$

この式の意味するところは、（混合ガウス分布に従って生成されると仮定した）データを１つ観測したときに、このデータはどのガウス分布から生成されたのか？を確率的に計算できる、ということ。

混合ガウス分布はkについてのシグマが増えた影響で、これまでのようにパラメータ推定の際に偏微分＝０　の方程式を解く方法が使えない。そこで、別の強力な推定方法にEMアルゴリズムがあるが、これは9章で。

2.4 指数型分布族, 2.4.2 共役事前分布

これまで二項分布の共役事前分布であるベータ分布、多項分布の共役事前分布であるディリクレ分布、…などいろいろな共役事前分布が紹介されてきた。では、共役事前分布は分布ごとに個別に分析して導出する必要があるのか？

...答えはNo。混合ガウス分布を除き、これまで扱ってきた確率分布はすべて以下の式の形で書くことが出来る指数型分布族 (exponential family) となっていた。

$$
p( \mathbf{x} | \mathbf{\eta} ) = h( \mathbf{x} ) g( \mathbf{\eta} )
\exp \{ \mathbf{ \eta }^T \mathbf{u}( \mathbf{x} ) \} \tag{2.194}
$$

※ベルヌーイ分布、多項分布、ガウス分布において$h(x), u(x), g(\eta)$の具体的な形を求める例がテキストにあるが略

指数型分布族である任意の分布はテキストの式(2.229)のように共役事前分布が導出できる。

2.4.3 無情報事前分布

（この節の話があまり分からない。以下、要注意）

ベイズ推論を行うことで、事後分布が得られる。すなわち、推定値がどの程度ばらつきそうか、という情報が得られる。一方で、ベイズ推論を行うには事前分布を決める必要がある。でも、事前分布に関する知見がない場合はどうすればいい？、というのがこの節のテーマ。

素朴な方法として、パラメータがとりうる範囲で一定値（一様分布）とすることを思いつく。でも、例えばガウス分布の平均の事前分布のように$ - \infty$から$\infty$まで取りうるパラメータの事前分布を一定値とする、つまりパラメータの事前分布$p(\lambda) = \text{const} $とすると、$ \int p(\lambda) d \lambda = \infty$となってしまう。

また、計算を簡単にするために変数変換を行うと、変数変換後の変数上では一定値でなくなってしまうという問題がある。

そこで、以下の2つの観点を満たすような事前分布はどのようなものかを考える。

平行移動不変性

これは事前分布の考える範囲を平行移動、例えば(0,10)の範囲で考えていたものを(-10,0)と移動しても、その範囲の確率が一定となる性質である。
（事前分布に関して知見がない前提なので、どの範囲を選んでも確率は一定となってほしいことに留意）

ガウス分布の平均$\mu$の無情報事前分布を考える場合、ガウス分布の平均の共役事前分布はガウス分布である。この事前分布のガウス分布の分散$\sigma^2_0$を$\infty$であると極限をとる（＝分布の「すそ」が無限に長く、極限に平らな形をしたガウス分布）といいらしい。

尺度不変性

これは事前分布の考える範囲をスケール変換、例えば(0,10)の範囲で考えていたものを(10, 100)や(100, 1000)としても、その範囲の確率が一定となる性質である。

平均$\mu$を考慮済みのガウス分布の$\sigma^2$の無情報事前分布を考える場合、計算の簡単化のために分散の代わりに精度$\lambda$に変数変換を行い、$\lambda$の事前分布として$a_0 = b_0 = 0$であるガンマ分布を使うといいらしい。

★分からないこと：平行移動不変性と尺度不変性の両方を満たすような無情報事前分布は与えられない・・？

2.5 ノンパラメトリック法

これまで扱ってきた確率分布は、例えば平均と共分散という２つのパラメータが決まると形が一意に決まるガウス分布であったり、ベクトル$\alpha$という１つのパラメータが決まると形が一意に決まるディリクレ分布、…のように少数のパラメータで形が決まるものを利用してきた。このアプローチはパラメトリックなアプローチと呼ばれる。

パラメトリックなアプローチの欠点は、これだと仮定した分布がデータを生成した真の分布を表現するのに不適切な場合はどれだけ頑張っても精度よくフィッティングできない。例えば、複数個の山があるような多峰性の分布は、単峰制のガウス分布では決して表現できない。

そこで、分布の形状にわずかな仮定しかしないノンパラメトリックな3つのアプローチに触れる。これから触れられる3アプローチは（ベイズ主義でなく）簡単な頻度主義的手法であり

• ヒストグラム密度推定法
• カーネル密度推定法（2.5.1節）
• 最近傍法（2.5.2節）

である。

ヒストグラム密度推定法

この手法は非常にシンプル。ヒストグラムの幅$\Delta$を適当に決め、その幅の間に入ったデータ数÷全データ数をその幅の確率とするというやり方（下図参照）。緑線が真の分布で、青棒のヒストグラムがヒストグラム密度推定法によって推定された確率密度分布。

ただし、デメリットが多いので、実務では1次元か2次元のデータの簡便な可視化に使うぐらいだろう。

＜デメリット＞

• 推定された密度（分布）の形は$\Delta$に大きく依存する（図参照）
• 推定した密度が区間の幅で不連続
• 次元数が増えると計算規模が指数的に増加（次元の呪い）

2.5.1 カーネル密度推定法

まずは、カーネル密度推定法および最近傍法の導出に使う確率密度の推定式である式(2.246)の導出を行おう。

未知の確率密度$p(x)$を推定したい。ある小さな領域Rに割り当てられる確率は

$$
P = \int_{R} p(x) dx \tag{2.242}
$$

で与えられる。
（1.2.1節「確率密度」で$x$が区間$(a,b)$にある確率は

$$
p ( x \in (a, b) ) = \int_a^b p(x) dx \tag{1.24}
$$

で与えられた。これの多次元へ拡張した自然な考え方）

ここで2.1節「二項分布」を思い出そう。表がでる確率$\mu$のコインを$N$回投げて$m$回表が出る確率は二項分布で

$$
\text{Bin}(m | N, \mu) = \frac{N!}{(N-m)! m!} \mu^m ( 1 - \mu)^{N-m} \tag{2.9}
$$

で表せた。
この考え方を使うと、領域Rに割り当てられる確率$P$な$N$個の観測データ集合のうち$K$個のデータがRに入る確率は

$$
\text{Bin}(K | N, P) = \frac{N!}{K! (N-K)!} P^K ( 1 - P)^{N-K} \tag{2.243}
$$

と書ける。

次に以下の期待値、分散の公式を紹介する。

• $a$を定数として$ \mathrm{E}[aX] = a \mathrm{E}[X] $
  （サイコロの目の期待値$\mathrm{E}[X]$は3.5、「サイコロの目の2倍」の期待値$\mathrm{E}[2X]$は
  「サイコロの目の期待値」の2倍の7と等しくなることをイメージすれば分かりやすい）
• $a$を定数として$ \text{var}[aX] = a^2 \text{var}[X] $
  （分散の計算式に2乗がでてくるので、$X$を$a$倍すると$a^2$がでてくるイメージ）

この公式と2.1節で導出した二項分布の期待値、分散の式を使うと、領域R内のデータ点$K$の期待値をNで割った$K/N$の期待値、および分散は以下のように書ける。

• $ \mathrm{E}[K/N] = P$
• $ \text{var}[K/N] = P (1-P) / N $

分散の式より$N$が大きくなるほど分散が小さくなり、$K/N$は$P$一点に固まることになる。なので、

$$
K \simeq NP \tag{2.244}
$$

と書ける。
次に領域Rが十分に小さく、確率密度$p(x)$がR内で一定とみなせるとすれば

$$
P \simeq p(x) V \tag{2.245}
$$

式(2.244), (2.245)より確率密度$p(x)$の推定値は

$$
p(x) = \frac{K}{NV} \tag{2.246}
$$

と考えることが出来る。

やっと準備ができたので、ここからカーネル密度推定法の議論を始めよう。
式(2.246)のVを固定し、Kをデータから推定する。すなわち、立方体やガウス関数で定義される領域Rの体積Vを固定とし、この領域内に含まれるデータ数Kを推定する。

領域内のデータ点を数えるために、以下の関数を定義する。
原点を中心とする辺の長さ１の単位立方体内にデータ$\mathbf{u}$があれば1、なければ0となるような関数。

\begin{align}

k( \mathbf{u} ) = 

\begin{cases}

1, & |u_i| \leq 1/2 \quad　( i= 1, ..., D ) \\
0, & else

\end{cases} \tag{2.247}
\end{align}

この式より、中心が$\mathbf{x}$で一辺が$h$の立方体内にデータ点$\mathbf{x} _ n$があれば１、なければ０となる関数は$k( (\mathbf{x} - \mathbf{x} _ n) / h)$と書ける。
よって、この立方体内の総データ数$K$は

$$
K = \sum_{n=1}^N k \left( \frac{ (\mathbf{x} - \mathbf{x} _ n) }{h} \right) \tag{2.248}
$$

となる。ここで、 $k( (\mathbf{x} - \mathbf{x} _ n) / h)$は$k( (\mathbf{x} _ n - \mathbf{x}) / h)$と考えても同じ（対称性）なので、２通りの見方ができる。

• 中心が$\mathbf{x}$で一辺が$h$の立方体内にデータ点$\mathbf{x}_n$があれば１、なければ０となる関数
  の総和
• 中心が$\mathbf{x} _n$にある立方体の総和（この文の意味が分かりにくいが、もう少し先まで読んで下の図をみると理解できる）

Vはいま立方体の体積である$h^D$で固定、Nはデータ数なので既知、Kは上式で推定したので確率密度の推定に必要な変数はすべて準備が出来た。式(2.246)のKに(2.248)を代入することで、確率密度は

$$
p( \mathbf{x} ) = \frac{1}{N} \sum_{n=1}^N \frac{1}{h^D}
k \left( \frac{\mathbf{x} - \mathbf{x} _ n}{h} \right) \tag{2.249}
$$

となった。この式は結局、下図のように観測した各データ点を中心とした立方体（図中の青線と橙線。辺の長さは分析者が設定するパラメータ）の足し合わせで定義される黒線を正規化したものを確率密度とすることを意味している。

ただ、式(2.249)の推定結果は、推定した確率密度が不連続というヒストグラム密度推定法の欠点は解決できていない。

そこで、より滑らかな関数として一般的に使われるガウス関数を使うと、立方体のときと同様の手順で以下の式が導出される。この式による推定結果は下図。

$$
p( \mathbf{x} ) = \frac{1}{N} \frac{1}{ (2 \pi h^2)^{D/2} } \exp
\left\{ - \frac{ || \mathbf{x}^2 - \mathbf{x} _n^2 || }{2 h^2} \right\}
\tag{2.250}
$$

ここでhはガウス分布の標準偏差で、分析者が設定する必要があるパラメータである。

各データ点を中心（平均）、分散は分析者が設定した値で固定としたガウス分布の足し合わせである下図黒線を正規化したものを確率密度とする、という解釈。

ガウス分布の標準偏差hを適切な値にしないと推定した確率密度がトゲトゲしすぎたり、滑らかになりすぎたりする。

2.5.2 最近傍法

テキストでは、特に式変形なく簡単な概念の紹介のみとなっている。なので、この記事でもさらっと流そう。

カーネル密度推定法では（ガウス分布の分散のような）幅を決めるパラメータ$h$がどのデータ点に対しても一定であることが問題であった。そのため、データ密度の高い領域では$h$が大きすぎると平滑化されすぎて、真の分布の構造がかき消されることがある。そのため、データ空間内の位置に応じて$h$を変えるようにしようとしているのが最近傍法の考え。

確率密度$p(x)$の推定値の式に戻る。

$$
p(x) = \frac{K}{NV} \tag{2.246}
$$

カーネル密度推定法ではVを固定しKの値を求めたが、最近傍法ではKを固定しVの値を求める。このために、$p(x)$を推定したい点$\mathbf{x}$を中心とした小球を考えて、ちょうどK個のデータ点を含むようになるまで半径を広げる。その方法で求めた確率密度分布は下図となり、Kが平滑化の度合いを決めていることが分かる。

1. 日常会話では「ノイズに対してロバスト」とかって言い方をよく聞く ↩