ぺ ん ぎ ん の 閃 き

閃き- blog

きらびやかに、美しく、痛烈に.

PAC学習と計算論的学習理論(Computatinonal Learning Theory)の文献まとめ

f:id:yumaloop:20190224225119p:plain


 PAC学習 (Probability Approximately Correct learning) とは、イギリスの理論計算機科学者 Leslie Valiant が1984年に以下の論文で初めて提唱した概念で、計算機科学の分野でそれまで研究されていた一般的な計算アルゴリズムの効率性・複雑性に対して、学習アルゴリズムに関しても同様の議論を行うために定義されたものです。

"A theory of the learnable"
Valiant, L. G. (1984). Communications of the ACM 1984 pp1134-1142.

 Valiant は、統計学の分野である「パターン認識(Pattern recognition)」や「決定理論(Decision theory)」の諸問題を、理論計算機科学の分野である計算複雑性理論(Computational Complexity Theory, CCT)と結びつけることで、学習アルゴリズムの持つ「複雑性」を定義しています。具体的に言うと、計算複雑性理論における「多項式時間で計算可能な解法が存在する問題(クラスP)」の定義と近しい形で、「実現可能な学習アルゴリズム」のクラスを定義しました。やがて、Valiantにより創始された学習アルゴリズムに関する一連の研究分野は計算論的学習理論(Computational Learning Theory, CLT)と呼ばれるようになり、今日に至ります。PAC学習の研究成果としては、例えば機械学習モデルの持つ汎化性能についての上界を与えるものなどがあり、近年ではソフトマージンSVMDeep Learningの汎化誤差に関して、計算論的学習理論の研究者から、たくさんの関連論文が出ています。

 2019年2月現在、機械学習の隆盛下においても、計算論的学習理論やPAC学習に関する日本語の解説はあまりありませんが、「PAC learning」でGoogle検索すると、Introduction, Tutorioal, Surveyといった名前のつく英語の解説文献がそこそこ見つかります。個人的には、以下の文献がわかりやすかったです。

"Overview of the Probably Approximately Correct (PAC) Learning Framework" (David Haussler, UC Santa Cruz, July 2010)

"PAC Learnability" (Karl Stratos, TTIC, September 2016)

 このポストでは、Valiantが導入した計算論的学習理論における重要な2つの概念"true error(generalization error)", "PAC learnability"の定義を紹介し、参考論文を備忘として留めておきます。

定義1(true error or generalization error)
 学習したい"概念"クラス  \mathcal{C} とそれに対応づけられた確率変数の集合  \mathcal{X} に対して、学習アルゴリズム  \mathcal{A}: X \in \mathcal{X} \mapsto c(X) \in \mathcal{C} を考える。ただし、  \mathcal{X}に対して定義される確率分布を  \mathcal{D} とおく。確率分布  \mathcal{D} に対して独立同分布となる  m コのサンプル  \left\{ x_1, \dots, x_m \right\} と、それに対応する"概念"  \left\{ c(x_1), \dots, c(x_m) \right\} のペアからなる訓練集合  S := { \bigl\{ x_i, c(x_i) \bigr\} }_{i=1}^{m} によって学習を行い仮説  h_{S}: \mathcal{X} \to C を立てる。このとき、仮説  h_{S}真の誤差(true error)あるいは汎化誤差(generalization error)を次式で定義する。

$$ error(h) := \underset{x ~ \sim~ \mathcal{D}}{Pr} \bigl( h_{S}(x) \neq c(x) \bigr) $$

定義2(PAC learnable)
 "概念"クラス  CPAC学習可能であるとは、以下の条件を満たす学習アルゴリズム  \mathcal{A} が存在し、かつ学習に必要な訓練集合の大きさ  m:=|S| の下限に対して  m \geq poly(1/\varepsilon, 1/\delta, dim(\mathcal{X}), size(c)) という関係を満たす多項式 poly(\cdot) が存在することである。

$$ \forall \varepsilon, \delta \in [0, 1], ~~ \forall m \in \mathbb{N}, ~~ \forall c \in C, ~~ \forall \mathcal{D} ~on~ \mathcal{X}, \\ \underset{S ~ \sim~ \mathcal{D}^{m} \times C^{m}}{Pr} \Bigl( \underset{x ~ \sim~ \mathcal{D}}{Pr} \bigl( h_{S}(x) \neq c \bigr) \leq \varepsilon \Bigr) \geq 1 - \delta $$


※ PAC学習では、ある訓練集合  S に対して、高い確率で「高い予測精度のモデル」を作ることができる学習アルゴリズム  \mathcal{A} を考える。"PAC Learnable"の定義は抽象的だが、確率分布によって定まる汎化誤差(true error)を念頭においており、PAC学習の分野の成果は、実際に使われている多くの確率モデル・学習アルゴリズムに対する誤差理論へ適用できる。

※ ある訓練集合  S に対して仮説と学習アルゴリズムの関係は  h_{S} \gets \mathcal{A}(S) で与えられることに注意。


 ところで、「学習アルゴリズム  A 」という表現は、計算機科学・計算理論の色が根強く反映されていると思います。すなわち、数理統計学(確率&解析学)のお気持ちとしては、学習は、「確率変数の集合間の写像(関数)を、有限集合から近似する計算操作」として解釈されますが、ここでは近似すべき関数の方が強調されているように感じます。

アラン・チューリングは「計算する」とは何か?を定義し「計算できること/できないこと」に厳密な線引きを与えましたが、当然「学習する」とは何か?にも定義が必要です。計算論的学習理論では、計算の枠組みを拡張することによって、学習を同様のスコープで捉えているように思います。


論文リスト

PAC Learningおよび計算論的理論に関する論文
PAC-Baysian Theory*1に関する論文

参考サイト

以下の記事は、PACを含む計算論的学習理論を概説しており、勉強になります。

qiita.com

*1:PAC LearningとBayse推定を組み合わせた理論

最尤推定量とワルド検定・スコア検定・尤度比検定

1. パラメータの尤もらしさに関する統計的仮説検定

 何らかのパラメトリックな確率モデル  f(x|\theta) とデータサンプル  X^{n} に対して定義される対数尤度関数  \ell(\theta | X^{n}) = \log f(X^{n}|\theta) に対して、「あるパラメータ  \theta最尤推定 \theta_{MLE} と一致しているかどうか」に真偽を与える統計的仮説検定を考えます。すなわち、仮説:

帰無仮説  H_0: \theta = \theta_{MLE}
・対立仮説  H_1: \theta \neq \theta_{MLE}


を設けて、帰無仮説  H_0: \theta = \theta_{MLE}有意水準  \alpha で棄却する条件を、データ(標本)に基づく統計量とその分布から求めます。

 このポストでは、上記の検定方式として、Wald検定Score検定尤度比検定の3つを紹介します。これらは、パラメータ空間  \Theta から  \mathbb{R}^{1} への写像である対数尤度関数  \ell(\theta | X^{n}) を中心に考えると、自然な導出であることがわかります。(下図参照)

f:id:yumaloop:20190206221533p:plain


1.1 ワルド検定(Wald test)

ワルド検定では、現在のパラメータ  \theta最尤推定 \theta_{MLE} とのを基に検定を行います。

f:id:yumaloop:20190206221834p:plain

最尤推定 \theta_{MLE} に対するワルド検定(Wald test)

最尤推定量の漸近正規性(Asymptotic normality):

$$ \frac{ \hat{\theta}_{MLE} - \theta }{ \sqrt{ \frac{1}{ I_{n}(\theta) } } } \underset{d}{\to} N \left( 0, 1 \right) $$

が成り立つとき、仮説:

帰無仮説  H_0: \theta = \theta_{MLE}
・対立仮説  H_1: \theta \neq \theta_{MLE}


における帰無仮説  H_0有意水準  \alpha の棄却域  R は、

$$ R = \left\{ x \in \mathcal{X} ~|~ \left| \frac{ \hat{\theta}_{MLE} - \theta }{ \sqrt{ \frac{1}{ I_{n}(\theta) } } } \right| \geq {Z}_{\frac{\alpha}{2}} \right\} $$

で与えられる。


1.2 スコア検定(Score test)

スコア検定では、現在のパラメータ  \theta の対数尤度関数における傾き(スコア関数)を基に検定を行います。

f:id:yumaloop:20190206221743p:plain

最尤推定 \theta_{MLE} に対するスコア検定(Score test)

スコア関数  S_n (\theta, X) とフィッシャ情報量  I_n(\theta) の定義

$$ \begin{eqnarray} S_n \left( \theta, X \right) & := & \frac{d}{d \theta} \log f_n (X | \theta) \\ I_n \left( \theta \right) & := & \mathbb{E} \left[ { S_n \left( \theta, X \right) }^{2} \right] \end{eqnarray} $$

より、スコア関数  S_n (\theta, X) の平均と分散が

$$ \begin{eqnarray} \mathbb{E} \left[ S_n \left( \theta, X \right) \right] & = & 0 \\ \mathbb{V} \left[ S_n \left( \theta, X \right) \right] & = & I_n (\theta) \end{eqnarray} $$

であることから、中心極限定理(CLT, Central Limit Theorem)を用いて、

$$ \frac{ S_n \left( \theta, X \right) }{ \sqrt{ I_{n}(\theta) } } \underset{d}{\to} N \left( 0, 1 \right) $$

が成り立つ。よって、仮説:

帰無仮説  H_0: \theta = \theta_{MLE}
・対立仮説  H_1: \theta \neq \theta_{MLE}


における帰無仮説  H_0有意水準  \alpha の棄却域  R は、

$$ R = \left\{ x \in \mathcal{X} ~|~ \left| \frac{ S_n \left( \theta, X \right) }{ \sqrt{ I_{n}(\theta) } } \right| \geq {Z}_{\frac{\alpha}{2}} \right\} $$

で与えられる。


1.3 尤度比検定(Likelihood ratio test)

尤度比検定では、現在のパラメータ  \theta の対数尤度と最尤推定 \theta_{MLE} との対数尤度のを基に検定を行います。

f:id:yumaloop:20190206221807p:plain

最尤推定 \theta_{MLE} に対する尤度比検定(Likelihood ratio test)

尤度比検定統計量  \lambda(X) を $$ \lambda(X) := \frac{ f_n (X | \theta)}{ f_n (X | \theta_{MLE}) } $$

と定義すると、仮説:

帰無仮説  H_0: \theta = \theta_{MLE}
・対立仮説  H_1: \theta \neq \theta_{MLE}


において、帰無仮説  H_0: \theta = \theta_{MLE} の下で、

$$ -2 \log \lambda(X) \underset{d}{\to} {\chi_1}^{2} $$

が成り立つことから、帰無仮説  H_0有意水準  \alpha の棄却域  R は、

$$ R = \left\{ x \in \mathcal{X} ~|~ -2 \log \lambda(X) \gt {\chi_{1, ~ \alpha}}^{2} \right\} $$

で与えられる。



2. KL-divergenceとFischer情報量の関係

 スコア関数  S_{n}\bigl( \theta, X^{n} \bigr) と、Fischer情報量  I_{n}\bigl( \theta \bigr) は以下のように定義されます。 $$ \begin{eqnarray} S_{n}\bigl( \theta, X^{n} \bigr) & := & \frac{d}{d \theta} \log f(X^{n}|\theta) \\ \\ I_{n}\bigl( \theta \bigr) & := & \mathbb{E}_{X^{n}~\sim~f(\cdot | \theta)} \left[ { \left( \frac{d}{d \theta} \log f(X^{n}|\theta) \right) }^{2} \right] \\ & = & \int_{- \infty}^{\infty} f(X^{n}|\theta) { \left( \frac{d}{d \theta} \log f(X^{n}|\theta) \right) }^{2} dx \end{eqnarray} $$

ここで、パラメータ  \theta の尤度  f(X^{n}|\theta) と、パラメータ  \theta の値を微小量  +h だけ変化させたときの尤度  f(X^{n}|\theta + h) との乖離度として、KL-divergenceを考えます。

$$ \begin{eqnarray} D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + h) \bigr) & = & \mathbb{E}_{X^{n}~\sim~f(\cdot | \theta)} \left[ \log \frac{f(X^{n}|\theta)}{f(X^{n}|\theta + h)} \right] \\ & = & \int_{- \infty}^{\infty} f(X^{n}|\theta) \log \frac{f(X^{n}|\theta)}{f(X^{n}|\theta + h)} dx \end{eqnarray} $$

さらに、  D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + h) \bigr) h の関数とみて、 h=0 のまわりでテイラー展開マクローリン展開)をすると、

$$ \begin{eqnarray} D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + h) \bigr) & = & D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + 0) \bigr) \\ & ~~~~~ + & \frac{1}{1!} \frac{d}{dh} D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + h) \bigr) {\large \mid}_{h=0} \cdot h^{1} \\ & ~~~~~ + & \frac{1}{2!} \frac{d^{2}}{dh^{2}} D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + h) \bigr) {\large \mid}_{h=0} \cdot h^{2} \\ & ~~~~~ + & O \left( h^{3} \right) \\ \\ & = & \mathbb{E}_{X^{n}} \left[ \log \frac{f(X^{n}|\theta)}{f(X^{n}|\theta)} \right] \\ & ~~~~~ + & \frac{1}{1!} \mathbb{E}_{X^{n}} \left[ \frac{d}{d \theta} \log f(X^{n}|\theta) \right] \cdot h^{1} \\ & ~~~~~ + & \frac{1}{2!} \mathbb{E}_{X^{n}} \left[ { \left( \frac{d}{d \theta} \log f(X^{n}|\theta) \right) }^{2} \right] \cdot h^{2} \\ & ~~~~~ + & O \left( h^{3} \right) \\ \\ & = & 0 \\ & ~~~~~ + & \mathbb{E}_{X^{n}} \left[ S_n \left( \theta, X^{n} \right) \right] \cdot h^{1} \\ & ~~~~~ + & \frac{1}{2} I_n \bigl( \theta \bigr) \cdot h^{2} \\ & ~~~~~ + & O \left( h^{3} \right) \\ \\ & = & \frac{1}{2} I_n \bigl( \theta \bigr) \cdot h^{2} + O \left( h^{3} \right) \end{eqnarray} $$

となります。つまり、 パラメータ  \theta の微小変化に対する KL-divergence: D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + h) \bigr) の局所二次近似は、パラメータ  \theta のFischer情報量と密接に関わります。

$$ \begin{eqnarray} D_{KL} \bigl( f( X^{n} | \theta) ~||~ f( X^{n} | \theta + h) \bigr) & = & \frac{1}{2} I_n \bigl( \theta \bigr) \cdot h^{2} + O \left( h^{3} \right) \\ & \simeq & \frac{1}{2} I_n \bigl( \theta \bigr) \cdot h^{2} \end{eqnarray} $$

以上の結果は、次の定理にまとめられます。

$$ \lim_{h \to 0} \frac{1}{h^{2}} D_{KL} \bigl( f(X^{n}|\theta) \mid\mid f(X^{n}|\theta+h) \bigr) = \frac{1}{2} I(\theta) $$

f:id:yumaloop:20190206235442p:plain

 最尤推定は、さまざまな確率モデルに対する最適化(学習)の基本であり、尤度関数(損失関数)の挙動を正確に把握することがとても重要です。Deep learningなど、高次元のパラメータ空間をもつ確率モデルの場合でも、二次の偏導関数(ヘッセ行列)によって、尤度関数の凸性・曲率が全て記述*1され、これはパラメータに対するFischer情報量行列と密に関連します。すなわち、Fischer情報量行列が正則行列でない*2場合、尤度関数はパラメータに対して非凸で複雑な関数となり、最尤推定量を探索することはより難しくなります。



3. 参考書籍

現代数理統計学 (創文社現代経済学選書)

現代数理統計学 (創文社現代経済学選書)

Statistics

Statistics

*1:線形空間では基底ベクトルに対する微分さえわかれば、あらゆる方向に対する勾配がその線形結合として表現可能です。

*2:すべての固有値が正である行列のことを、正則行列といいます。正則行列ではない行列を特異行列ということがあります。

ガウス過程と回帰モデル(線形~線形基底~ガウス過程)

f:id:yumaloop:20190203151221g:plain



※ このポストの内容は、PRML下巻6章とGPML2章を参考にしています。

www.microsoft.com www.gaussianprocess.org

1. ガウス過程(GP, Gaussian Process)

1.1 ガウス分布の共役性

 任意の確率モデル  p(\cdot|w) に対して  w の事前分布  p(w)としてガウス分布を仮定した場合、ベイズの定理によって求められる  w の事後分布  p(w|X^{n})ガウス分布に従います。つまり、ガウス分布に従うパラメータ  w について、任意の確率モデル  p(x|w) とデータサンプル  X^{n} から、事後分布  p(w|X^{n}) を求める操作はガウス分布の中で閉じています。一般に、「事前分布としてある確率分布を仮定すると、事後分布も同じ確率分布に従う」という性質を共役性といい、ガウス分布は共役性をもちます。

$$ \begin{eqnarray} (Posterior~of~w) & = & \frac{(likelihood ~of~ X^{n})\times(Prior~of~w)}{(Evidence~of~the~model)} \\ \\ p(w|X^{n}) & = & \frac{p(X^{n}|w)p(w)}{\int p(X^{n}|w)p(w) ~dw} \end{eqnarray} $$ $$ p(w) \sim N(\mu, s) ~ \Rightarrow ~ p(w|X^{n}) \sim N(\mu', s') $$

1.2 ガウス分布の切断と周辺化

 それぞれがガウス分布に従う「互いに独立でない」確率変数列  Y = \left\{ Y_1, \dots Y_N \right\} を考えます。このとき、確率変数列  Y はNコの確率変数  \left\{ Y_i \right\}_{i=1 \sim N} のペアであるため、確率変数列  Y そのものを  \mathbb{R}^{N} 次元の確率変数として扱うことができます。同様に、確率変数列  Y における「任意の  m コの部分列」は  \mathbb{R}^{m} 次元の確率変数として扱えます。

 平均が  \mu で、共分散行列が  \Sigma のN次元ガウス分布に従う確率変数  X

 
  X \in \mathbb{R}^{N}, ~~~ X \sim N_{N} \bigl( \mu, \Sigma \bigr)


に対して、 X_a, X_b a+b=N)による分割を考えます。

 X = \left[ \begin{array}{c} X_a \\ X_b \end{array} \right] \sim  N_{a+b} \left( \left[ \begin{array}{c} \mu_a \\\ \mu_b \end{array} \right], ~ \left[ \begin{array}{c} \Sigma_{aa} \\ \Sigma_{ba} \end{array} \begin{array}{c} \Sigma_{ab} \\ \Sigma_{bb} \end{array} \right] \right) \\\
 
  \mu = \left[ \begin{array}{c} \mu_a \\\ \mu_b \end{array} \right], ~~~
  \Sigma = \left[ \begin{array}{c} \Sigma_{aa} \\ \Sigma_{ba} \end{array} \begin{array}{c} \Sigma_{ab} \\ \Sigma_{bb} \end{array} \right]


このとき、 X_a, X_b は以下のようなガウス分布に従います。


  X_{a} \sim N_a \left( \mu_a, \Sigma_{aa} \right) \\
  X_{b} \sim N_b \left( \mu_b, \Sigma_{bb} \right)


さらに、「 X の共分散行列  \Sigma は半正定値対称行列である」ことから、次式が成り立ちます。


  \Sigma_{ab} = \Sigma_{ba}^{-1}


ここで、「 X_a, X_b は独立ではない」ことに注意すると、条件つき確率  p(X_a | X_b) は次のように与えられ、条件つき平均  \mu_{a|b} と条件つき共分散行列  \Sigma_{a|b} からなるガウス分布に従うことがわかります。


  X_{a} | X_{b} \sim N \left( \mu_{a|b}, \Sigma_{a|b} \right) \\

\begin{eqnarray}
  \left\{
    \begin{array}{c}
      \mu_{a|b} & = & \mu_a + \Sigma_{ba}{\Sigma_{aa}}^{-1} (X_b - \mu_b) \\
      \Sigma_{a|b} & = & \Sigma_{aa} - \Sigma_{ab} {\Sigma_{bb}}^{-1} \Sigma_{ba}
    \end{array}
  \right.
\end{eqnarray}


 X_a を手元にあるデータサンプル、  X_b を新しいデータサンプルと捉えると、これらを結合した  X_{a+b} に対してガウス過程を導入することができます。

1.3 ガウス過程の定義

 「任意の組が同時ガウス分布(joint Gaussian distribution)に従う」ような確率変数の集合をガウス過程(Gaussian Process, GP)と定義します。ガウス分布ガウス過程を比較すると以下のようになります。

 
  {\bf X} \sim N \left( \mu, \Sigma \right) as: \\\
 
  \begin{eqnarray}
    {\bf \mu} & := &  \mathbb{E} \left[ {\bf X} \right] \\
    \Sigma & := & \mathbb{C o v} \left[ X_i, X_j) \right] = \mathbb{E} \bigl[ X_i - \mu_i \bigr] \mathbb{E} \bigl[ X_j - \mu_j \bigr]
  \end{eqnarray}


  • ガウス過程(Gaussian Process, GP):任意のガウス過程  \mathcal{GP}() は、平均関数  m({\bf x}) と共分散関数  k({\bf x}, {\bf x}') によって完全に記述される。

  f({\bf x}) \sim \mathcal{GP} \left( m({\bf x}), ~ k({\bf x}, {\bf x}') \right) as: \\\

\begin{eqnarray}
  m({\bf x}) & := & \mathbb{E} \left[ f({\bf x}) \right] \\\
  k({\bf x}, {\bf x}') & := & \mathbb{C o v} \left[ f({\bf x}), f({\bf x}') \right] = \mathbb{E} \bigl[ f({\bf x} - m({\bf x}) \bigr] \mathbb{E} \bigl[ f({\bf x}' - m({\bf x}') \bigr]
\end{eqnarray}




2. さまざまな回帰モデルの比較

2.1 回帰問題とは?

 まず、回帰(Regression)に一般的な定義を与えてみます。

 確率変数  X \in \mathbb{R}^{d}, Y \in \mathbb{R} の実現値からなるNコのデータサンプル:

$$ \mathcal{D} := \left\{ \left( \begin{array}{c} {\bf x}_1 \\ y_1 \end{array} \right), \dots, \left( \begin{array}{c} {\bf x}_N \\ y_N \end{array} \right) \right\} $$

を用いて、 X の値から  Y の値を推定(予測)するモデル:

$$ f: X \in \mathbb{R}^{d} \to Y \in \mathbb{R} $$

を作りたい。

ここで、モデル  f による推定値  f({\bf x}_i) と真値  y_i との推定誤差  \varepsilon_i を次のように定義し、

$$ \varepsilon_i := y_i - f({\bf x}_i) ~~~ (i=1, \dots, N) $$

さらに、各  \varepsilon_i ~({\small i=1, \dots, N}) はそれぞれ独立に同じパラメータをもつガウス分布に従うと仮定すると、

$$ y_i = f({\bf x}_i) + {\varepsilon}_i ~~~ (i=1, \dots N) \\ {\varepsilon}_i \sim N \left( 0, \frac{1}{\beta} \right) $$

という等式が得られます。これをまとめると、

$$ p( y_i | {\bf x}_i, f) = N \left( f({\bf x}_i), \frac{1}{\beta} \right) $$

となり、回帰問題の目的は、「ばらつきを持つ確率分布  p({y | x_i, f}) を定めること」としてベイズ的に解釈することも可能です。すなわち、次式(ベイズの定理)で、  p( f | y_i, {\bf x}_i ) はデータサンプル  \mathcal{D} への「当てはまりの良さ(=尤度)」を表しますが、同時に、  p( f | y_i, {\bf x}_i ) によって、データサンプル  \mathcal{D} から  f の事後確率と予測分布が決定されます。

$$ p( f | y_i, {\bf x}_i ) = \frac{ p( y_i | {\bf x}_i, f ) p(f) } {p( y_i | {\bf x}_i )} $$

 一般化に、回帰問題では、データサンプル  \mathcal{D} = \left\{ {\bf X}, {\bf y} \right\} とモデル  f に対して

$$ p( f | {\bf y}, {\bf X} ) = \frac{ p( {\bf y} | {\bf X}, f ) p(f) } {p( {\bf y} | {\bf X} )}, ~~ i.e. ~~ p( f | \mathcal{D} ) = \frac{ p( \mathcal{D} | f ) p(f) } {p( \mathcal{D} )} $$

という構造が現れます。


2.2 線形回帰モデル

2.2.1 線形回帰モデルの準備

 もっとも自然な回帰モデルとして「線形回帰モデル」を考えます。モデルパラメータとして、重みベクトル  {\bf w} \in \mathbb{R}^{d} を想定します。

$$ y_i = f({\bf x}_i) + \varepsilon_i, \\ f({\bf x}_i) = {\bf x}_i^{\rm T} {\bf w} = \left( \begin{array}{c} x_{i1}, \dots, x_{id} \end{array} \right) \cdot \left( \begin{array}{c} w_{1} \\ \vdots \\ w_{d} \end{array} \right), ~~~ \varepsilon_i \sim N \left( 0, \frac{1}{\beta} \right) $$

これらの関係式をまとめると、各  {\bf x}_i に対する  y_i の推定値は、次式のように「ばらつきを持った確率モデル」として表現されます。

$$ p( y_i | x_i, {\bf w} ) = N \left( y_i | f({\bf x}_i), \frac{1}{\beta}{\bf I} \right) = N \left( {\bf x}_i^{\rm T} {\bf w}, \frac{1}{\beta} \right) $$

さらに、個々データではなくデータセット  \mathcal{D} 全体に対して、回帰モデル  f はどのように作用するのかをまとめて考えます。計画行列(design matrix)  {\bf X} \in \mathbb{R}^{d ~ \times ~ N} と、ベクトル  {\bf y} を以下のように定義します。

$$ {\bf X} := \left( {\bf x}_1, {\bf x}_2, \dots, {\bf x}_N \right)= \left( \begin{array}{cccc} x_{11} & x_{21} & \dots & x_{N1}\\ x_{12} & x_{22} & \dots & x_{N2}\\ \vdots & \vdots & \dots & \vdots \\ x_{1N} & x_{2N} & \dots & x_{NN} \end{array} \right), ~~~ {\bf y} := \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N \end{array} \right), ~~~ {\boldsymbol \varepsilon} := \left( \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{array} \right) $$

すると、「線形回帰モデル」は以下のように記述され、

$$ {\bf y} = f({\bf X}) + {\boldsymbol \varepsilon}, \\ f({\bf X}) = {\bf X}^{\rm T}{\bf w}, ~~~ \boldsymbol \varepsilon \sim N \left( {\bf 0}, \frac{1}{\beta}{\bf I}\right) $$

これらの関係式をまとめると、 データセット  \mathcal{D} における  {\bf y} の推定値は、次式のように「ばらつきを持った確率モデル」として表現されます。

$$ p( {\bf y} | {\bf X}, {\bf w} ) = N \left( {\bf y} | f({\bf X}), \frac{1}{\beta}{\bf I} \right) = N \left( {\bf X}^{\rm T} {\bf w}, \frac{1}{\beta}{\bf I} \right) $$

2.2.2 線形回帰モデルのベイズ的解釈

 線形回帰モデルにおいて、 {\bf y} を推定する確率モデル  p( {\bf y} | {\bf X}, {\bf w} ) の性能はパラメータ  {\bf w} の値によって左右されます。確率モデル  p( {\bf y} | {\bf X}, {\bf w} )ベイズの枠組み:

$$ p({\bf w} | {\bf y}, {\bf X} ) = \frac{ p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) } {p({\bf y} | {\bf X} )} = \frac{ p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) } { \int p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) ~ d{\bf w}} $$

によって拡張し、パラメータ  {\bf w} の事後分布  p({\bf w} | {\bf y}, {\bf X} ) を導出することで、N+1 番目の未知のデータサンプル  \left( \begin{array}{c} {\bf x}^{*} \\ y^{*} \end{array} \right) に対する「ばらつきを持った確率モデル」である予測分布を導出しましょう。

(1) {\bf w} の事前分布

$$ p({\bf w}) = N \left( {\bf w} ~ | ~ {\bf 0}, ~ \frac{1}{\alpha} {\bf I} \right) $$

(2)尤度  p({\bf y} | X, {\bf w})

$$ \begin{eqnarray} p({\bf y} | X, {\bf w}) & = &\prod_{i=1}^{N} p(y_i | {\bf x}_i, {\bf w}) \\ & = & \prod_{i=1}^{N} \frac{1}{ \sqrt{ 2 \pi \beta^{-1}}} \exp \left( - \frac{{\left( y_i - {{\bf x}_i}^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & \frac{1}{ {\left( 2 \pi \beta^{-1} \right)}^{n/2}} \exp \left( - \frac{{\left( {\bf y} - X^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & N \left( {\bf y} ~ | ~ X^{\rm T}{\bf w}, ~ \frac{1}{\beta}{\bf I} \right) \end{eqnarray} $$

(3) {\bf w} の事後分布

$$ \begin{eqnarray} p({\bf w} | {\bf X}, {\bf y}) & \propto & \exp \left( - \frac{1}{2} { ( {\bf y} - {\bf X}^{\rm T}{\bf w} ) }^{\rm T} { \left( \frac{1}{\beta} {\bf I} \right) }^{-1} ( {\bf y} - {\bf X}^{\rm T}{\bf w} ) \right) \exp \left( - \frac{1}{2} {\bf w}^{\rm T} { \left( \frac{1}{\alpha} {\bf I} \right) }^{-1} {\bf w} \right) \\ & \propto & \exp \left( { ( {\bf w} - {\hat{\bf w}}_{MAP} ) }^{\rm T} \left( \beta {\bf X}{\bf X}^{\rm T} + \alpha {\bf I} \right) ( {\bf w} - {\hat{\bf w}}_{MAP} ) \right) \\ & = & N \left( {\bf w} ~ | ~ {\hat{\bf w}}_{MAP}, ~ {\bf A}^{-1} \right) \end{eqnarray} $$

 ただし、

$$ \begin{eqnarray} {\hat{\bf w}}_{MAP} & := & \beta {\bf A }^{-1} {\bf X}{\bf y} \\ {\bf A} & := & \beta {\bf X}{\bf X}^{\rm T} + \alpha {\bf I} \end{eqnarray} $$

(4) y^{*} の予測分布

$$ \begin{eqnarray} p^{*}(y^{*} | {\bf x}^{*}, {\bf X}, {\bf y}) & = & \int p(y^{*} | {\bf x}^{*}, {\bf w}) p( {\bf w} | X, {\bf y}) ~ d{\bf w} \\ & = & N \left( y^{*} ~|~ \beta {{\bf x}^{*}}^{\rm T} { \bf A }^{-1} {\bf X} {\bf y}, ~ { {\bf x}^{*}}^{\rm T} { \bf A }^{-1} {\bf x}^{*} \right) \end{eqnarray} $$

2.3 線形基底モデル

2.3.1 線形基底モデルの準備

 線形モデルでは、入力ベクトル  {\bf x} \in \mathbb{R}^d のあるデータ空間(data space) \subset \mathbb{R}^d の中で、出力  y \in \mathbb{R} を推定します。しかし、線形回帰モデルではデータ空間を構成する基底ベクトルの線形結合によって  y を推定するため、モデルの表現力が制約されてします。そこで、入力ベクトル  {\bf x} を基底関数(basis function)  \phi: \mathbb{R}^d \to \mathbb{R}^M によって特徴量  \phi({\bf x}) \in \mathbb{R}^M へと変換し、高次元の特徴空間  \subset \mathbb{R}^{M}(feature space)から  y \in \mathbb{R} の値を推定することを考えます。ここで、基底関数  \phi による変換は線形でも非線形でもよいとします。

$$ {\bf x}_i \in \mathbb{R}^d \to \phi({\bf x}_i) \in \mathbb{R}^M \underset{\large \bf w}{\longrightarrow} f({\bf x}_i) \in \mathbb{R} \underset{ {\large \varepsilon_i} } {\longleftrightarrow} y_i \in \mathbb{R}, ~~~ (i=1, \dots, N) \\ $$

基底関数  \phi: \mathbb{R}^d \to \mathbb{R}^M は、データを表現する基底を変える(空間を変える)関数と考えることができます。また、一般に、 d \lt M であることから、 \phi はデータサンプルを高次元空間に写像する(mapping)関数であることがわかります。

  •  {\bf x}_i のデータ空間  \subset \mathbb{R}^d へのmapping: $$ {\bf x}_i = \left( x_{i1}, x_{i2}, \dots, x_{id} \right) $$
  •  {\bf x}_i の特徴空間  \subset \mathbb{R}^M へのmapping: $$ \phi({\bf x}_i) = \left( \phi_{1}({\bf x}_i), \phi_{2}({\bf x}_i), \dots, \phi_{M}({\bf x}_i) \right) $$

 「線形基底モデル」を考えます。モデルパラメータとして、重みベクトル  {\bf w} \in \mathbb{R}^{M} を想定します。

$$ y_i = f({\bf x}_i) + \varepsilon_i \\ f({\bf x}_i) = { \phi ( {\bf x}_i ) }^{\rm T} {\bf w} = \left( \begin{array}{c} \phi_{1}({\bf x}_i), \dots, \phi_{M}({\bf x}_i) \end{array} \right) \cdot \left( \begin{array}{c} w_{1} \\ \vdots \\ w_{M} \end{array} \right), ~~~ \varepsilon_i \sim N \left( 0, \frac{1}{\beta} \right) $$

これらの関係式をまとめると、各  {\bf x}_i に対する  y_i の推定値は、次式のように「ばらつきを持った確率モデル」として表現されます。

$$ p( y_i | x_i, {\bf w} ) = N \left( y_i | f( {\bf x}_i ), \frac{1}{\beta} \right) = N \left( y_i | {\phi ( {\bf x}_i ) }^{\rm T} {\bf w}, \frac{1}{\beta} \right) $$

さらに、個々データではなくデータセット  \mathcal{D} 全体に対して、回帰モデル  f はどのように作用するのかをまとめて考えます。計画行列(design matrix)  \Phi({\bf X}) \in \mathbb{R}^{M ~ \times ~ N} と、ベクトル  {\bf y}, {\boldsymbol \varepsilon} を以下のように定義します。

$$ \Phi({\bf X}) := \left( \phi({\bf x}_1), \phi({\bf x}_2), \dots, \phi({\bf x}_N) \right)= \left( \begin{array}{cccc} \phi_{1}({\bf x}_1) & \phi_{1}({\bf x}_2) & \dots & \phi_{1}({\bf x}_N) \\ \phi_{2}({\bf x}_1) & \phi_{2}({\bf x}_2) & \dots & \phi_{2}({\bf x}_N) \\ \vdots & \vdots & \dots & \vdots \\ \phi_{M}({\bf x}_1) & \phi_{M}({\bf x}_2) & \dots & \phi_{M}({\bf x}_N) \end{array} \right), ~~~ {\bf y} := \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N \end{array} \right), ~~~ {\boldsymbol \varepsilon} := \left( \begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{array} \right) $$

すると、「線形基底回帰モデル」は以下のように記述され、

$$ {\bf y} = f({\bf X}) + {\boldsymbol \varepsilon}, \\ f({\bf X}) = {\Phi({\bf X})}^{\rm T}{\bf w}, ~~~ \boldsymbol \varepsilon \sim N \left( {\bf 0}, \frac{1}{\beta}{\bf I}\right) $$

これらの関係式をまとめると、 データセット  \mathcal{D} における  {\bf y} の推定値は、次式のように「ばらつきを持った確率モデル」として表現されます。

$$ p( {\bf y} | {\bf X}, {\bf w} ) = N \left( {\bf y} ~|~ f({\bf X}), \frac{1}{\beta}{\bf I} \right) = N \left( {\bf y} ~|~ {\Phi({\bf X})}^{\rm T}{\bf w}, \frac{1}{\beta}{\bf I} \right) $$

2.3.2 線形基底モデルのベイズ的解釈

 確率モデル  p( {\bf y} | {\bf X}, {\bf w} )ベイズの枠組み:

$$ p({\bf w} | {\bf y}, {\bf X} ) = \frac{ p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) } {p({\bf y} | {\bf X} )} = \frac{ p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) } { \int p({\bf y} | {\bf X}, {\bf w}) p({\bf w}) ~ d{\bf w}} $$

によって拡張し、パラメータ  {\bf w} の事後分布  p({\bf w} | {\bf y}, {\bf X} ) を導出することで、N+1 番目の未知のデータサンプル  \left( \begin{array}{c} {\bf x}^{*} \\ y^{*} \end{array} \right) に対する「ばらつきを持った確率モデル」である予測分布を導出しましょう。

(1) {\bf w} の事前分布

$$ p({\bf w}) = N \left( {\bf w} ~ | ~ {\bf 0}, ~ \frac{1}{\alpha} {\bf I} \right) $$

(2)尤度  p({\bf y} | X, {\bf w})

$$ \begin{eqnarray} p({\bf y} | \Phi({\bf X}), {\bf w}) & = &\prod_{i=1}^{N} p(y_i | \phi({\bf x}_i), {\bf w}) \\ & = & \prod_{i=1}^{N} \frac{1}{ \sqrt{ 2 \pi \beta^{-1}}} \exp \left( - \frac{{\left( y_i - {\phi({\bf x}_i)}^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & \frac{1}{ {\left( 2 \pi \beta^{-1} \right)}^{n/2}} \exp \left( - \frac{{\left( {\bf y} - {\Phi({\bf X})}^{\rm T}{\bf w} \right)}^{2}}{2 \beta^{-1}} \right) \\ & = & N \left( {\bf y} ~ | ~ {\Phi({\bf X})}^{\rm T}{\bf w}, ~ \frac{1}{\beta}{\bf I} \right) \end{eqnarray} $$

(3) {\bf w} の事後分布  p({\bf w} | \Phi({\bf X}), {\bf y})

$$ \begin{eqnarray} p({\bf w} | \Phi({\bf X}), {\bf y}) & \propto & \exp \left( - \frac{1}{2} { ( {\bf y} - {\Phi({\bf X})}^{\rm T}{\bf w} ) }^{\rm T} { \left( \frac{1}{\beta} {\bf I} \right) }^{-1} ( {\bf y} - {\Phi({\bf X})}^{\rm T}{\bf w} ) \right) \exp \left( - \frac{1}{2} {\bf w}^{\rm T} { \left( \frac{1}{\alpha} {\bf I} \right) }^{-1} {\bf w} \right) \\ & \propto & \exp \left( { ( {\bf w} - {\hat{\bf w}}_{MAP} ) }^{\rm T} \left( \beta {\Phi({\bf X})}{\Phi({\bf X})}^{\rm T} + \alpha {\bf I} \right) ( {\bf w} - {\hat{\bf w}}_{MAP} ) \right) \\ & = & N \left( {\bf w} ~ | ~ {\hat{\bf w}}_{MAP}, ~ { A }^{-1} \right) \end{eqnarray} $$

 ただし、

$$ \begin{eqnarray} {\hat{\bf w}}_{MAP} & := & \beta { \bf A }^{-1} {\bf X}{\bf y} \\ A & := & \beta {\Phi({\bf X})}{\Phi({\bf X})}^{\rm T} + \alpha {\bf I} \end{eqnarray} $$

(4) y^{*} の予測分布

$$ \begin{eqnarray} p^{*}(y^{*} | {\bf x}^{*}, \Phi({\bf X}), {\bf y}) & = & \int p(y^{*} | {\bf x}^{*}, {\bf w}) p( {\bf w} | \Phi({\bf X}), {\bf y}) ~ d{\bf w} \\ & = & N \left( y^{*} ~|~ \beta {{\bf x}^{*}}^{\rm T} { \bf A}^{-1} \Phi({\bf X}) {\bf y}, ~ { {\bf x}^{*}}^{\rm T} { \bf A}^{-1} {\bf x}^{*} \right) \end{eqnarray} $$

2.4 ガウス過程による線形基底モデルの表現

2.4.1 ガウス過程による回帰モデルの準備

 ガウス過程による回帰では、基底関数  \phi() ではなく、回帰モデル  f({\bf X}) の確率分布  p(f({\bf X})) = p( {\Phi({\bf X})}^{\rm T}{\bf w} ) の共分散行列  K_{XX} を事前に定めます。  K_{XX} をグラム行列(Gram Matrix)といいます。

$$ {\bf K}_{XX} = \left( \begin{array}{cccc} k({\bf x}_1, {\bf x}_1) & k({\bf x}_1, {\bf x}_2) & \dots & k({\bf x}_1, {\bf x}_N) \\ k({\bf x}_2, {\bf x}_1) & k({\bf x}_2, {\bf x}_2) & \dots & k({\bf x}_2, {\bf x}_N) \\ \vdots & \vdots & \ddots & \vdots \\ k({\bf x}_N, {\bf x}_1) & k({\bf x}_N, {\bf x}_2) & \dots & k({\bf x}_N, {\bf x}_N) \\ \end{array} \right) $$ ただし、 $$ \begin{eqnarray} k({\bf x}_i, {\bf x}_j) & = & C o v \left[ f( {\bf x}_i ), f( {\bf x}_j ) \right] \\ & = & { \phi( {\bf x}_i )}^{\rm T} \left( \frac{1}{\alpha} {\bf I} \right) \phi( {\bf x}_j ) \\ & = & \frac{1}{\alpha} \phi_{1}( {\bf x}_i ) \phi_{1}( {\bf x}_j ) + \cdots + \frac{1}{\alpha} \phi_{M}( {\bf x}_i ) \phi_{M}( {\bf x}_j ) \end{eqnarray} $$

 定義から、グラム行列  {\bf K}_{XX} f({\bf x}_i) f({\bf x}_j) の分散からなる行列(=共分散行列)です。すなわち、グラム行列  {\bf K}_{XX} は、基底関数(basis function)  \phi() と計画行列(design matrix)  \Phi({\bf X}) および重みベクトル  {\bf w} の事前分布  p(w) \sim N \left( {\bf 0}, \frac{1}{\alpha}{\bf I} \right) によっても記述されます。

$$ \begin{eqnarray} {\bf K}_{XX} & = & \left( \begin{array}{cccc} C o v \left[ f({\bf x}_1), f({\bf x}_1) \right] & C o v \left[ f({\bf x}_1), f({\bf x}_2) \right] & \dots & C o v \left[ f({\bf x}_1), f({\bf x}_N) \right] \\ C o v \left[ f({\bf x}_2), f({\bf x}_1) \right] & C o v \left[ f({\bf x}_2), f({\bf x}_2) \right] & \dots & C o v \left[ f({\bf x}_2), f({\bf x}_N) \right] \\ \vdots & \vdots & \ddots & \vdots \\ C o v \left[ f({\bf x}_N), f({\bf x}_1) \right] & C o v \left[ f({\bf x}_N), f({\bf x}_2) \right] & \dots & C o v \left[ f({\bf x}_N), f({\bf x}_N) \right] \\ \end{array} \right) \\ \\ & = & {\Phi({\bf X})}^{\rm T} \left( \frac{1}{\alpha}{\bf I} \right) \Phi({\bf X}) \\ & = & \frac{1}{\alpha} \Phi({\bf X}) { \Phi({\bf X}) }^{\rm T} \end{eqnarray} $$

よって、線形基底モデルを「①基底関数  \phi()、計画行列  \Phi({\bf X}) によって表現する(weight space view)」「②グラム行列  {\bf K}_{XX} によって表現する(function space view)」ことが可能です。

$$ \begin{eqnarray} p \bigl( f({\bf X}) \bigr) & = & N \left( f({\bf X}) ~|~ {\bf 0}, \frac{1}{\alpha}\Phi({\bf X}) {\Phi({\bf X})}^{\rm T} \right) \\ & = & N \Bigl( f({\bf X}) ~|~ {\bf 0}, {\bf K}_{XX} \Bigr) \\ \\ p \bigl( {\bf y} \bigr) & = & N \left( {\bf y} ~|~ {\bf 0}, \frac{1}{\alpha}\Phi({\bf X}) {\Phi({\bf X})}^{\rm T} + \frac{1}{\beta} {\bf I} \right) \\ & = & N \Bigl( {\bf y} ~|~ {\bf 0}, {\bf K}_{XX} + \frac{1}{\beta} {\bf I} \Bigr) \\ \end{eqnarray} $$

2.4.2 ガウス過程による回帰モデルによる予測

 ここで、N+1 番目の未知のデータ  {\bf x}^{*} に対する  y^{*} の予測値  f( {\bf x}^{*} ) を導出しましょう。

$$ y^{*} = f( {\bf x}^{*} ) + {\varepsilon}^{*} \\ f( {\bf x}^{*} ) = { \phi( {\bf x}^{*} ) }^{\rm T} {\bf w}, ~~~ {\varepsilon}^{*} \sim N \left( 0, \frac{1}{\beta} \right) $$

2.3.2 の結果から、 f ( {\bf x}^{*} ) の予測分布は以下のように記述されます。

$$ f ( {\bf x}^{*} ) \sim N \Bigl( \mathbb{E} \left[ f ( {\bf x}^{*} ) \right] , \mathbb{V} \left[ f ( {\bf x}^{*} ) \right] \Bigr) $$ $$ \begin{eqnarray} \left\{ \begin{array}{c} \mathbb{E} \left[ f ( {\bf x}^{*} ) \right] & = & { {\bf k}^{*} }^{\rm T} { \left( {\bf K}_{XX} + \frac{1}{\beta} {\bf I} \right) }^{-1} {\bf y} \\ \mathbb{V} \left[ f ( {\bf x}^{*} ) \right] & = & k ( {\bf x}^{*}, {\bf x}^{*} ) - { {\bf k}^{*} }^{\rm T} { \left( {\bf K}_{XX} + \frac{1}{\beta} {\bf I} \right) }^{-1} {\bf k}^{*} \end{array} \right. \end{eqnarray} $$

未知のデータ  {\bf x}^{*} ガウス過程  \mathcal{GP} \left( {\bf 0}, {\bf K}_{XX} \right) に加えると、

$$ \begin{eqnarray} {{\bf K}_{XX}}^{*} & = & \left( \begin{array}{cccc} k( {\bf x}_1, {\bf x}_1 ) & \dots & k( {\bf x}_1, {\bf x}_N ) & k( {\bf x}_1, {\bf x}^{*} ) \\ \vdots & \ddots & \vdots & \vdots \\ k( {\bf x}_N, {\bf x}_1 ) & \dots & k( {\bf x}_N, {\bf x}_N ) & k( {\bf x}_N, {\bf x}^{*} ) \\ k( {\bf x}^{*}, {\bf x}_1 ) & \dots & k( {\bf x}^{*}, {\bf x}_N ) & k( {\bf x}^{*}, {\bf x}^{*} ) \\ \end{array} \right) \\ \\ & = & \left( \begin{array}{cc} {\bf K}_{XX} & {\bf k}^{*} \\ { {\bf k}^{*} }^{\rm T} & k( {\bf x}^{*}, {\bf x}^{*} ) \end{array} \right) \end{eqnarray} $$

ただし、 {\bf k}^{*} は以下のように記述されるベクトルです。 $$ \begin{eqnarray} {\bf k}^{*} & := & \left( \begin{array}{c} k( {\bf x}_1, {\bf x}^{*} ) \\ \vdots \\ k( {\bf x}_N, {\bf x}^{*} ) \end{array} \right) \\ \\ & = & \left( \begin{array}{c} C o v \left[ f( {\bf x}_1 ), f( {\bf x}^{*} ) \right] \\ \vdots \\ C o v \left[ f( {\bf x}_N ), f( {\bf x}^{*} ) \right] \end{array} \right) \\ \\ & = & \left( \begin{array}{c} { \phi( {\bf x}_1 )}^{\rm T} \left( \frac{1}{\alpha} {\bf I} \right) \phi( {\bf x}^{*} ) \\ \vdots \\ { \phi( {\bf x}_N )}^{\rm T} \left( \frac{1}{\alpha} {\bf I} \right) \phi( {\bf x}^{*} ) \end{array} \right) \\ \\ & = & \left( \begin{array}{c} \frac{1}{\alpha} \phi_{1}( {\bf x}_1 ) \phi_{1}( {\bf x}^{*} ) + \cdots + \frac{1}{\alpha} \phi_{M}( {\bf x}_1 ) \phi_{M}( {\bf x}^{*} ) \\ \vdots \\ \frac{1}{\alpha} \phi_{1}( {\bf x}_N ) \phi_{1}( {\bf x}^{*} ) + \cdots + \frac{1}{\alpha} \phi_{M}( {\bf x}_N ) \phi_{M}( {\bf x}^{*} ) \end{array} \right) \\ \\ & = & { \phi( {\bf x}^{*} ) }^{\rm T} \frac{1}{\alpha} {\bf I} ~ \Phi( {\bf X} ) \end{eqnarray} $$



3. ガウス過程回帰をPythonで実装

 実際に、ある関数からデータサンプルを作成して、ガウス過程による回帰モデル  f によってこれを近似します。今回は、以下のような平均関数  m(x) と共分散関数  k(x, x') を用いたガウス過程  \mathcal{GP} と、それに従う回帰モデル  f を考えます。ここで、共分散関数としてRBFカーネル(Radial basis function kernel, Gaussian kernel)を採用しています。

$$ f(x) \sim \mathcal{GP} \left( m(x), k(x, x') \right) $$

 
\begin{eqnarray}
  \left\{
    \begin{array}{c}
      m(x) & = & 0 \\
      k(x, x') & = & \exp \left( - \frac{ { || x - x' || }^{2} }{2 {\sigma}^{2} } \right) = \exp \left( - \gamma { || x - x' || }^{2}  \right)
    \end{array}
  \right.
\end{eqnarray}


この回帰モデルを用いて、関数:

$$ y = \sin(2 \pi x) + \cos(1.7 \times 2 \pi x) ~~~ (0 \leq x \leq 1) $$

を、生成されるデータサンプルから予測してみましょう。なお実装は、こちらのQuitta記事を参考にさせてもらいました。

3.1 RBFカーネルのためのクラス

class RBFKernel(object):
    
    def __init__(self, sigma):
        self.sigma = sigma
        
    def get_sigma(self):
        return np.copy(self.sigma)

    def __call__(self, x1, x2):
        return np.exp( (-1. / (self.sigma ** 2) ) *  (x1 - x2) ** 2)

    def derivatives(self, x1, x2):
        dif_sigma = np.exp( (-1. / (self.sigma ** 2) ) *  (x1 - x2) ** 2) * ( (x1 - x2) ** 2 ) / ( self.sigma ** 3)
        return dif_sigma

    def update_sigma(self, update):
        self.sigma += update

3.2 ガウス過程回帰を行うためのクラス

class GaussianProcessRegression(object):

    def __init__(self, kernel, x, y, beta=1.):
        self.kernel = kernel
        self.x = x
        self.y = y
        self.beta = beta
        
        # グラム行列を計算
        self.fit_kernel()

    def fit_kernel(self):
        ''' グラム行列を計算 '''
        
        # グラム行列 K
        self.gram = self.kernel(*np.meshgrid(self.x, self.x))
        # 共分散行列 K + 1/beta I
        self.covariance = self.gram + np.identity(len(self.x)) / self.beta
        # 精度行列 (共分散行列の逆行列)
        self.precision = np.linalg.inv(self.covariance)
    
    def predict_y(self, x_new):
        ''' y_new の予測値と分散を出力 '''
        
        # k*
        k = np.array([self.kernel(xi, x_new) for xi in self.x])

        # 平均と分散
        y_mean = k.dot(self.precision.dot(self.y))
        y_var = self.kernel(x_new, x_new) - k.dot(self.precision.dot(k))
        
        return y_mean, y_var

3.3 グラフの描画(100フレームのアニメーション)を作成する

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def func(x):
    ''' 回帰モデルによって近似を行う関数 '''
    return np.sin(2 * np.pi * x) + np.cos(1.7 * 2 * np.pi * x)


def create_data_sample(func, low=0, high=1., n=10, std=1.):
    ''' n 組のデータサンプル (x, y) を作成する '''
    x = np.random.uniform(low, high, n)
    y = func(x) + np.random.normal(scale=std, size=n)
    return x, y
    
def _update(frame, x_data, y_data):
    ''' フレームごとにグラフを描画する '''
        
    # 現在のグラフを消去する
    plt.cla()
    
    # データサンプルのframe数分の取得
    x = x_data[:frame]
    y = y_data[:frame]
    
    # 回帰モデルをつくる
    kernel = RBFKernel(sigma=0.5)
    regression = GaussianProcessRegression(kernel, x, y, beta=100)
    
    # 描画用にxのデータを生成
    x_test = np.linspace(0, 1, 100)
    
    # yの予測値の平均と分散を導出
    y_test=[]
    y_var=[]
    for x_new in x_test:
        y_new_mean, y_new_var = regression.predict_y(x_new)
        y_test.append(y_new_mean)
        y_var.append(y_new_var)
        
    # list を numpy.ndarrayに変換
    y_test = np.array(y_test)
    y_var = np.array(y_var)
    
    # グラフを描画
    plt.scatter(x, y, alpha=0.4, color="blue", label="data sample")
    plt.plot(x_test, func(x_test), color="blue", linestyle="dotted", label="actual function")
    plt.plot(x_test, y_test, color="red", label="predicted function")
    plt.fill_between(x_test, y_test - 3 * np.sqrt(y_var), y_test + 3 * np.sqrt(y_var), color="pink", alpha=0.5, label="95% confidence interval (2 * std)")
    plt.legend(loc="lower left")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.xlim([-0, 1])
    plt.ylim([-2.1, 2.1])
    plt.title("Gaussian Process Regression \n (RBF kernel, sigma=0.5, sample size = " + str(frame) + ")")
    
def main():
    # 描画領域
    fig = plt.figure(figsize=(8, 6))

    # データサンプルを作成
    x_data, y_data = create_toy_data(func, low=0, high=1, n=100, std=0.1)

    params = {
        'fig': fig,
        'func': _update,  # グラフを更新する関数
        'fargs': (x_data, y_data),  # 関数の引数 (フレーム番号を除く)
        'interval': 1,  # 更新間隔 (ミリ秒)
        'frames': np.arange(1, 101, 1),  # フレーム番号を生成するイテレータ, 1~100
        'repeat': False,  # 繰り返さない
    }

    # アニメーションを作成する
    ani = animation.FuncAnimation(**params)

    # アニメーションを保存する
    ani.save("./gpr_rbf.gif")
    # アニメーションを表示する
    # plt.show()

3.4 結果:アニメーションの描画

   サンプルサイズ(データ点の数)が増えると、徐々に回帰モデルで予測した関数(predicted function)が真の関数(actual function)に近づいていることがわかります。

f:id:yumaloop:20190203151221g:plain



4. 参考書籍

Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning series)

Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning series)

  • 作者: Carl Edward Rasmussen,Christopher K. I. Williams,Francis Bach
  • 出版社/メーカー: The MIT Press
  • 発売日: 2005/11/23
  • メディア: ハードカバー
  • 購入: 1人 クリック: 3回
  • この商品を含むブログ (2件) を見る

Pattern Recognition and Machine Learning (Information Science and Statistics)

Pattern Recognition and Machine Learning (Information Science and Statistics)

カーネル多変量解析―非線形データ解析の新しい展開 (シリーズ確率と情報の科学)

カーネル多変量解析―非線形データ解析の新しい展開 (シリーズ確率と情報の科学)

指数分布とワイブル分布をPythonでプロットしてみる


 生存時間解析など、応用範囲の広い指数分布についてまとめます。指数型分布族の仲間としては、ワイブル分布・ガンマ分布の他にも、ポアソン分布、レイリー分布、ラプラス分布(両側指数分布)などがあります。

1. 前提

1.1 確率分布の定義

 指数分布確率密度関数  f と累積分布関数  F は、次のように定義されます。

  • 指数分布

$$ X \sim Exp(\lambda), ~~~ \lambda \geq 0 $$

$$ \begin{eqnarray} f(x | \lambda) & = & \lambda ~ exp \left( - \lambda x \right) ~~~ (x \geq 0) \\ F(x | \lambda) & = & 1 - exp \left( - \lambda x \right) ~~~ (x \geq 0) \\ \end{eqnarray} $$ $$ \begin{eqnarray} \mathbb{E}\left[ X \right] & = & \frac{1}{\lambda} \\ \mathbb{V}\left[ X \right] & = & \frac{1}{\lambda^{2}} \end{eqnarray} $$

 また、ワイブル分布は、次のように定義されます。なお、 \Gamma(~) は、ガンマ関数*1(階乗を一般化したもの)です。

  • ワイブル分布

$$ X \sim Weibull(\alpha, \beta), ~~~ \alpha \geq 0, ~~~ \beta \geq 0 $$

$$ \begin{eqnarray} f(x | \alpha, \beta) & = & \frac{\alpha}{\beta} { \left( \frac{x}{\beta} \right) }^{\alpha - 1} \left\{ - { \left( \frac{x}{\beta} \right) }^{\alpha} \right\} ~~~ (x \geq 0) \\ F(x | \alpha, \beta) & = & 1 - exp\left\{ - {\left( \frac{x}{\beta} \right) }^{a} \right\} ~~~ (x \geq 0) \\ \end{eqnarray} $$ $$ \begin{eqnarray} \mathbb{E}\left[ X \right] & = & \beta ~ \Gamma \left( 1 + \frac{1}{\alpha} \right) \\ \mathbb{V}\left[ X \right] & = & \beta^{\large \frac{2}{\alpha}} \left\{ \Gamma(1 + \frac{2}{\alpha}) - {\left( \Gamma(1 + \frac{1}{\alpha}) \right)}^{2} \right\} \end{eqnarray} $$


1. 2 指数分布とワイブル分布の関係

 ワイブル分布のパラメータをそれぞれ  \alpha=1,~ \beta=\frac{1}{\lambda}とすると指数分布に一致します。つまり、ワイブル分布は、指数分布の一般に拡張したものになります。

$$ Exp(\lambda) = Weibull(1, \frac{1}{\lambda}) $$

ちなみに、以下で定義されるガンマ分布のパラメータをそれぞれ  \alpha=1,~ \beta=\frac{1}{\lambda}とすると指数分布に一致します。

$$ Exp(\lambda) = Weibull(1, \frac{1}{\lambda}) = Gamma(1, \frac{1}{\lambda}) $$

  • ガンマ分布

$$ X \sim Gamma(\alpha, \beta) $$

$$ \begin{eqnarray} f(x | \alpha, \beta) & = & \frac{1}{\Gamma(\alpha)} \frac{1}{\beta} {\left( \frac{x}{\beta} \right) }^{\alpha - 1}~~~ (x \gt 0) \\ \end{eqnarray} $$ $$ \begin{eqnarray} \mathbb{E}\left[ X \right] & = & \alpha \beta \\ \mathbb{V}\left[ X \right] & = & \alpha {\beta}^{2} \end{eqnarray} $$


2. Pythonによる実装

2.1 指数分布

python で分布を定義して、グラフをプロットしてみます。

import math
import random
import numpy as np
import matplotlib.pyplot as plt

# 指数分布

class Exp():
    def __init__(self):
        pass
    
    # 密度関数f
    def f(self, x, la=1):
        return la * math.exp(-1 * la * x)
    
    # 分布関数F
    def F(self, x, la=1):
        return 1 - math.exp(-1 * la * x)
    
    # Fの逆関数(samplingに使う)
    def invF(self, y, la=1):
        return -1 * (1 / la) * math.log(1 - y)
    
    # サンプリング
    def sampling(self, num, la):
        sample=[]
        for i in range(num):
            rand = random.random()
            sample.append(self.invF(rand, la))
        return np.array(sample)
    
    # 散布図
    def scatter_pdf(self, ax, data, la, alpha=0.4):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, la))

        ax.scatter(data, pdf, alpha=alpha)
        ax.set_title("Exponential (lambda={:.2f})".format(la))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)
        
    # グラフ
    def plot_pdf(self, ax, data, la):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, la))

        ax.plot(data, pdf)
        ax.set_title("Exponential (lambda={:.2f})".format(la))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)

実際にグラフを描画してみます。パラメータ  \lambda の値を変えて、分布の形状をみます。

# グラフをmatplotlibで描画してみる

data = np.linspace(0, 15, 100)

fig, ax = plt.subplots(3, 3, figsize=(18, 18))

Exp().plot_pdf(ax[0, 0], data, 0.2)
Exp().plot_pdf(ax[0, 1], data, 0.4)
Exp().plot_pdf(ax[0, 2], data, 0.6)
Exp().plot_pdf(ax[1, 0], data, 0.8)
Exp().plot_pdf(ax[1, 1], data, 1)
Exp().plot_pdf(ax[1, 2], data, 1.2)
Exp().plot_pdf(ax[2, 0], data, 1.4)
Exp().plot_pdf(ax[2, 1], data, 1.6)
Exp().plot_pdf(ax[2, 2], data, 1.8)

plt.show()
f:id:yumaloop:20190114000052p:plain:w1000

 プロットをみると、指数分布のパラメータ  \lambda は、尺度(ばらつき)を表していることがわかります。分布は  \lambda が大きくなるにつれ急峻になり、ばらつきが少なくなっていきます。平均について  \frac{1}{\lambda} \to 0 ~ (\lambda \to \infty)、分散について  \frac{1}{\lambda^{2}} \to 0 ~ (\lambda \to \infty) となることから、「指数分布はパラメータ  \lambda を上げると、平均は0に接近し、分散は限りなく小さくなる(0になる)」ということがわかります。

2.2 ワイブル分布

python で分布を定義して、グラフをプロットしてみます。

import math
import random
import numpy as np
import matplotlib.pyplot as plt

# ワイブル分布

class Weibull():
    def __init__(self):
        pass
        
    # 密度関数f
    def f(self, x, a, b):
        return (a/b) * ( (x/b)**(a - 1) ) * math.exp(-1 * (x/b)**a )
    
    # 分布関数F
    def F(self, x, a, b):
        return 1 - math.exp(-1 * (x/b)**a )
    
    # Fの逆関数
    def invF(self, y, a, b):
        return b * ( -1 * math.log(1 - y) )**(1/a)
    
    # サンプリング
    def sampling(self, num, a, b):
        sample=[]
        for i in range(num):
            rand = random.random()
            sample.append(self.invF(rand, a, b))
        return np.array(sample)
    
    # 散布図
    def scatter_pdf(self, ax, data, a, b, alpha=0.4):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, a, b))

        ax.scatter(data, pdf, alpha=alpha)
        ax.set_title("Weibull (alpha={:.2f}, beta={:.2f})".format(a, b))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)
        
    # グラフ
    def plot_pdf(self, ax, data, a, b):
        pdf=[]
        for d in data:
            pdf.append(self.f(d, a, b))

        ax.plot(data, pdf)
        ax.set_title("Weibull (alpha={:.2f}, beta={:.2f})".format(a, b))
        ax.set_xlabel("value")
        ax.set_ylabel("probability density")
        ax.set_ylim(0, 1)

実際にグラフを描画してみます。パラメータ  \alpha, ~\beta の値を変えて、分布の形状をみます。

# グラフをmatplotlibで描画してみる

data = np.linspace(0, 15, 100)

fig, ax = plt.subplots(3, 3, figsize=(18, 18))

Weibull().plot_pdf(ax[0, 0], data, 1, 1)
Weibull().plot_pdf(ax[0, 1], data, 1, 2)
Weibull().plot_pdf(ax[0, 2], data, 1, 3)
Weibull().plot_pdf(ax[1, 0], data, 2, 1)
Weibull().plot_pdf(ax[1, 1], data, 2, 2)
Weibull().plot_pdf(ax[1, 2], data, 2, 3)
Weibull().plot_pdf(ax[2, 0], data, 3, 1)
Weibull().plot_pdf(ax[2, 1], data, 3, 2)
Weibull().plot_pdf(ax[2, 2], data, 3, 3)

plt.show()
f:id:yumaloop:20190114000114p:plain:w1000

 プロットをみると、ワイブル分布のパラメータ  \alpha は分布の形状を表し、パラメータ  \beta は、尺度(ばらつき)を表していることがわかります。実際、分布の形状は  \alpha によって支配されており、ワイブル分布は  \alpha=1 のとき指数分布、  \alpha=2 のときレイリー分布と一致します。分布は  \beta が大きくなるにつれ急峻になり、ばらつきが少なくなっていきます。

*1:ガンマ関数は、$$ \Gamma(n):=\int_{0}^{\infty} y^{n - 1} e^{-y} dy $$と定義され、 n自然数ならば \Gamma(n)=(n-1)! が成り立つ。

Kullback-Leibler Divergenceについてまとめる


1. KL-divergenceとは?

 統計学情報理論をはじめとした広い分野で、KL-divergence*1はたびたび登場します。KL-divergenceは、「尤度比(尤もらしさを比較する尺度)を log 変換し(乗算操作を線形結合に直す、凸関数だから最適化との相性良い)、期待値(確率密度の重み付きの積分)をとったもの」として定義されます。

1.1 定義

 一般に、確率分布  P, ~ Q確率密度関数  p(x), ~ q(x) をもつとき、KL-divergence(Kullback-Leibler divergence)は、以下のように定義されます。

$$ D_{KL}( Q \mid\mid P ) := \int q(x) \log \frac{q(x)}{p(x)} ~dx \tag{1} $$


1.2 基本的な性質

 KL-divergenceは以下のような性質をもちます。

  • 非負性)非負の値域をもつ。

$$ 0 \leq D_{KL}( Q \mid\mid P ) \leq \infty \tag{2.1} $$

  • 完備性)値が 0 のとき、2つの確率分布 P と Q は同等。

$$ D_{KL}( Q \mid\mid P ) = 0 ~~ \Leftrightarrow ~~ P = Q \tag{2.2} $$

  • 非対称性)値は P と Q に対して対称性を持たない。

$$ D_{KL}( Q \mid\mid P ) \neq D_{KL}( P \mid\mid Q ) \tag{2.3} $$

  • 絶対連続性)値が発散しない限り、Q は P に関して絶対連続である。

$$ D_{KL}( Q \mid\mid P ) \lt \infty ~~ \Rightarrow ~~ P \gg Q \tag{2.4} $$


 例として、2つの正規分布間でKL-divergenceを計算*2すると、以下のような値になります。分布のズレが大きくなるほど、値が増加していることがわかります。

f:id:yumaloop:20190109162028p:plain:w500


1.3 KL-divergenceは距離なのか?

 KL-divergenceは、確率や情報を考える際にとても重要な量であるため、分野や文脈によって様々な名称で呼ばれます。

"KL-divergence"
"KL-metrics"
"KL-information"
"Information divergence"
"Information gain"
"Relative entropy"

 KL-divergenceはつねに非負の値をとるため、これは確率分布 P と Q が存在する空間における距離(metrics)を示していると解釈することができます。しかし、KL-divergenceは次の距離の公理のうち「非負性」と「完備性」しか満たさないため、厳密には 距離(metrics) ではありません。

距離  d(~) の公理

  1. 非負性     d(x, ~ y) \geq 0
  2. 完備性     d(x, ~ y) = 0 ~~ \Leftrightarrow ~~ x = y
  3. 対称性     d(x, ~ y) = d(y, ~ x)
  4. 三角不等式   d(x, ~ y) + d(y, ~ z) \geq d(x, ~ z)

 例えば、ユークリッド距離、二乗距離、マハラビノス距離、ハミング距離などは、この4つの条件をすべて満たすため、明らかに距離(metrics)として扱うことができます。一方で、KL-divergenceは距離(metrics)ではなくdivergenceです。"divergence"とは、距離の公理における4つの条件のうち「非負性」「対称性」のみを採用したものであり、"距離(metrics)"の拡張概念です。距離の公理による制約を減らし、抽象度の高い議論をしようというわけです。

 なお、このdivergenceという語は発散と訳されることが一般的で、例えば物理学ではベクトル作用素 div として登場します。divergenceのイメージに該当する日本語はありませんが、相違度、分離度、逸脱度、乖離度などが使われているようです。

 実例として、2つの正規分布  N(0, 1) (青)と  N(1, 2) (赤)の間のKL-divergenceを測ってみましょう。左が青からみた赤とのKL距離、右が赤からみた青とのKL距離となっており値が異なることがわかります。なお「等分散正規分布間のKL-divergence」という特殊な条件下では、対称性が成り立ちます。

f:id:yumaloop:20190115092949p:plain


2. 諸量との関係

2.1 KL-divergenceと相互情報量

 エントロピー  H(X)結合エントロピー  H(X,Y) 、条件つきエントロピー  H(X|Y)相互情報量  MI(X,Y) は、確率密度  Pr(~) を用いて以下のように定義されます。*3

$$ \begin{eqnarray} H(X) & := & - \int Pr(x) \log Pr(x) ~dx \tag{3.1} \\ H(X,Y) & := & - \int Pr(x,y) \log Pr(x,y) ~dy~dx \tag{3.2} \\ H(X|Y) & := & - \int Pr(x,y) \log Pr(x|y) ~dx~dy \tag{3.3} \\ MI(X,Y) & := & \int \int Pr(x,y) \log \frac{Pr(x,y)}{Pr(x)Pr(y)} ~dxdy \tag{3.4} \end{eqnarray} $$

2つの確率変数 X, Y について、相互情報量 MI(: Mutual Information)によって相互の関係性が明示されます。

f:id:yumaloop:20190107163949p:plain:w400

$$ \begin{eqnarray} MI(X,Y) & = & H(X) - H(X|Y) \tag{3.5.1} \\ & = & H(Y) - H(Y|X) \tag{3.5.2} \\ & = & H(X) + H(Y) - H(X,Y) \tag{3.5.3} \end{eqnarray} $$

 ここで、KL-divergenceと相互情報量 MI には以下の関係が成り立ちます。つまり、相互情報量 MI(X, Y) は、「確率変数 X と Y が独立でない場合の同時分布 P(X, Y)」と「独立である場合の同時分布 P(X)P(Y)」 の離れ具合(平均的な乖離度合い)を示しているのだと解釈できます。

$$ \begin{eqnarray} MI(X, Y) & = & D_{KL} \bigl( Pr(x, y) \mid\mid Pr(x)Pr(y) \bigr) \tag{3.6.1} \\ & = & \mathbb{E}_{Y} \bigl[ D_{KL} \bigl( Pr(x|y) \mid\mid Pr(x) \bigr) \bigr] \tag{3.6.2} \\ & = & \mathbb{E}_{X} \bigl[ D_{KL} \bigl( Pr(y|x) \mid\mid Pr(y) \bigr) \bigr] \tag{3.6.3} \\ \end{eqnarray} $$

(例)3.6.2の式展開

$$ \begin{eqnarray} MI(X,Y) & = & \int \int Pr(x,y) \log \frac{Pr(x,y)}{Pr(x)Pr(y)} ~dxdy \\ & = & \int \int Pr(x|y)Pr(y) \log \frac{Pr(x|y)Pr(y)}{Pr(x)Pr(y)} ~dxdy \\ & = & \int Pr(y) \int Pr(x|y) \log \frac{Pr(x|y)}{Pr(x)} ~dx~dy \\ & = & \int Pr(y) \cdot D_{KL} \bigl( Pr(x|y) \mid\mid Pr(x) \bigr) ~dy \\ & = & \mathbb{E}_{Y} \bigl[ D_{KL} \bigl( Pr(x|y) \mid\mid Pr(x) \bigr) \bigr] \end{eqnarray} $$


2.2 KL-divergenceと対数尤度比

 ベイズ推定や統計モデリングでは、しばしば「真の分布  q を、 p(x|\hat{w})(パラメータの推定値  \hat{w} と確率モデル  p(x|w) の組み合わせ)推定する」というモチベーションが起こります。そのため、2つの分布のズレを測りたいとき、あるいは推定誤差を損失関数やリスク関数に組み込んでパラメータに対する最適化問題をときたいときに、KL-divergenceを使います。

 KL-divergenceは最尤法をベースに、尤度比検定(Likelihood ratio test)、ベイズ因子(Bayse factor)、赤池情報量基準(AIC, Akaike's Information Critation)などのモデル選択手法*4と深い関わりを持ちます。

  •  真の分布  P に対する推定分布  p のKL-divergence: D_{KL}(Q \mid\mid P) は次式のように、「  q p の対数尤度比に対して、真の分布  q で平均操作をおこなった値」です。

$$ \left( 対数尤度比 \right) = \log \frac{q(x)}{p(x)} \tag{4} $$

$$ \begin{eqnarray} D_{KL}( Q \mid\mid P ) & := & \int q(x) \log \frac{q(x)}{p(x)} ~dx \\ & = & \mathbb{E}_{X} \left[ \log \frac{q(x)}{p(x)} \right] \left(: 平均対数尤度比 \right) \tag{5.1} \\ \end{eqnarray} $$

  モデル比較・選択の評価値としてKL-divergenceを用いる場合、次式のように「KL-divergenceを最小にすること」と「平均対数尤度を最大にする」ことは等価になります。

$$ \begin{eqnarray} D_{KL}( Q \mid\mid P ) & = & \mathbb{E}_{X} \bigl[ \log q(x) \bigr] - \mathbb{E}_{X} \bigl[ \log p(x) \bigr] \tag{5.2} \\ & \propto & - \mathbb{E}_{X} \bigl[ \log p(x) \bigr] \left(: - 平均対数尤度 \right) \tag{5.3} \end{eqnarray} $$


  • パラメトリックな確率モデル  f(x|\theta) (例えば線形回帰モデル)を考えた場合、何らかの損失関数  L に対して最適なパラメータ  \theta_{0} と、その推定値  \hat{\theta} に対して、確率モデルの誤差は以下のようにKL-divergenceで表現されます。(ただし、 \ell( \cdot| x ) を対数尤度関数とする。)

$$ \left( 対数尤度比 \right) = \log \frac{f(x|\theta_{0})}{f(x|\hat{\theta})} \tag{6} $$

$$ \begin{eqnarray} \hat{\theta} & :=& \underset{\theta \in \Theta}{\rm argmin} ~ L(\theta) \tag{7} \\ \\ D_{KL}( \theta_{0} \mid\mid \hat{\theta} ) & := & \int f(x|\theta_{0}) \log \frac{f(x|\theta_{0})}{f(x|\hat{\theta})} ~dx \\ & = & \mathbb{E}_{X} \left[ \log \frac{ f(x|\theta_{0}) }{ f(x|\hat{\theta}) } \right] \tag{8.1} \\ & = & \mathbb{E}_{X} \bigl[ \ell( \theta_{0}|x ) \bigr] - \mathbb{E}_{X} \bigl[ \ell( \hat{\theta}|x ) \bigr] \tag{8.2} \end{eqnarray} $$


2.3 KL-divergenceとFisher情報量

 確率モデル  f(\cdot | \theta) を仮定したとき、パラメータ  \thetaFisher情報量  I(\theta) は以下のように定義されます。(ただし、 \ell( \cdot | x ) を対数尤度関数とする。)

$$ \begin{eqnarray} I(\theta) & := & \mathbb{E}_{X} \left[ { \left\{ \frac{d}{dx}\ell(\theta | x) \right\} }^{2} \right] \tag{9.1} \\ & = & \mathbb{E}_{X} \left[ { \left\{ \frac{d}{dx}\log f(x|\theta) \right\} }^{2} \right] \tag{9.2} \end{eqnarray} $$    確率モデル  f(\cdot | \theta) を仮定したとき、KL-divergenceとFisher情報量  I(\theta) の間には、次のような関係*5が成り立ちます。

$$ \begin{eqnarray} \lim_{h \to 0} \frac{1}{h^{2}} D_{KL} \bigl( f(x|\theta) \mid\mid f(x|\theta+h) \bigr) = \frac{1}{2} I(\theta) \tag{10} \end{eqnarray} $$

 この式は、 \theta の近傍でのKL-divergence: D_{KL} \bigl( f(x|\theta) \mid\mid f(x|\theta+h) \bigr) は、 \theta における確率モデル  f(\cdot|\theta) のFisher情報量  I(\theta) に比例することを示しています。つまり、Fisher情報量  I(\theta) は確率モデル  f(\cdot | \theta) \theta の近傍でもつ局所的な情報量を表していることがわかります。

f:id:yumaloop:20190107190651p:plain:w600


3. 参考書籍

数理統計学 (数学シリーズ)

数理統計学 (数学シリーズ)

Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing)

Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing)

情報量規準 (シリーズ・予測と発見の科学)

情報量規準 (シリーズ・予測と発見の科学)



*1:一般化するとf-divergenceというクラスが定義できます。汎関数微分などの関数解析の知識が必要。。

*2:pythonのライブラリ関数scipy.stats.entropyを使いました。

*3:熱力学的なエントロピーはボルツマン発祥ですが、シャノンの情報量に関する歴史的背景については以下で言及されています。ハートレーナイキスト→シャノンという流れがあるみたい。http://www.ieice.org/jpn/books/kaishikiji/200112/200112-9.html

*4:情報量規準についての文献→ https://www.ism.ac.jp/editsec/toukei/pdf/47-2-375.pdf

*5:この式は、KL-divergenceをテイラー展開して二次形式で近似すると導出されます。$$ \ell(\theta + h) - \ell(\theta) = {\ell}^{'}(\theta)h + \frac{1}{2} {\ell}^{''}(\theta) h^{2} + O(h^{3}) $$

最尤推定・MAP推定・ベイズ推定を比較する


1. 推定のモチベーション

 確率変数  X の真の分布  q(X) に対して、

$$ \begin{cases} 確率モデル  : p(X|w) \\ 事前分布   : \varphi(w) \end{cases} $$

を仮定したとき、

$$ データサンプル   : X^{n} := \left\{ X_1, \dots ,X_n \right\} $$

を用いて、 q(X) を推定したい。



2. 最尤推定(ML)

2.1 パラメータ  w の対数尤度  \ell(w|\cdot)

ベイズの定理: $$ \begin{eqnarray} p(w | X^{n}) = \frac{ p(X^{n}|w) \varphi(w) }{ \int_{W} p(X^{n}|w) \varphi(w) dw} \end{eqnarray} $$

において、以下で定義される対数尤度  \ell(w|\cdot) は、手元にあるサンプルデータ  X^{n}に対する 、確率モデル  p(x|w) の尤もらしさを示します。対数尤度  \ell(w|\cdot) の示す値は、サンプルデータ  X^{n} がもつ"ばらつき"が加味されておらず、手元にあるデータの値に大きく依存するということに注意してください。

$$ \begin{eqnarray} \ell(w | \cdot) & := & \log p(\cdot | w) \\ \\ \ell(w | X^{n}) & := & \log p(X^{n} | w) \\ & = & \log \prod_{i=1}^{n} p(X_i | w) \end{eqnarray} $$

2.2 パラメータの推定(最尤推定量)

 最尤推定では、対数尤度  \ell が最大となるようなパラメータ  \hat{w}_{ML} を選び、これを確率モデル  p(X|w) におけるパラメータ  w の推定量とします。推定量  \hat{w}_{ML}最尤推定と呼びます。

$$ \begin{eqnarray} \hat{w}_{ML} := \underset{w \in W}{\rm argmax} ~ \ell(w|X^{n}) \end{eqnarray} $$

以上から、最尤推定では、確率変数  X の真の分布  q(x) は、サンプルデータ  X^{n} から求めた最尤推定  \hat{w}_{ML}確率モデル  p(X|w) によって次式のように推定されます。

$$ \begin{eqnarray} q(X) \approx p( X | \hat{w}_{ML} ) \end{eqnarray} $$



3. 事後確率最大化推定(MAP)

3.1 パラメータ  w の事後確率  p(w|\cdot)

ベイズの定理:

$$ \begin{eqnarray} p(w | X^{n}) = \frac{ p(X^{n}|w) \varphi(w) }{ \int_{W} p(X^{n}|w) \varphi(w) dw} \end{eqnarray} $$

において、以下で定義される事後確率  p(w|\cdot) は、手元にあるサンプルデータ  X^{n}に対する 、確率モデル  p(X|w)事前分布  \varphi(w) の尤もらしさを示します。事後確率  p(w|\cdot) の示す値は、サンプルデータ  X^{n} がもつ"ばらつき"が加味されておらず、その値に大きく依存するということに注意してください。

$$ \begin{eqnarray} p(w | \cdot) & := & \frac{1}{Z} \log p(\cdot | w)\varphi(w) \\ \\ p(w | X^{n}) & := & \frac{1}{Z_n} \log p(X^{n} | w)\varphi(w) \\ & = & \frac{1}{Z_n} \log \prod_{i=1}^{n} p(X_i | w)\varphi(w) \end{eqnarray} $$


  • ただし、 Z,~Z_n は、それぞれベイズの定理における右辺の分子に対して  w積分した値(正規化定数)を表す。

3.2 パラメータの推定(MAP推定量

 MAP推定では、事後確率が最大となるようなパラメータ  \hat{w}_{MAP} を選び、これを確率モデル  p(X|w) におけるパラメータ  w の推定量とします。推定量  \hat{w}_{MAP}MAP推定量と呼びます。

$$ \begin{eqnarray} \hat{w}_{MAP} := \underset{w \in W}{\rm argmax} ~ p(w|X^{n}) \end{eqnarray} $$

 以上から、MAP推定では、確率変数  X の真の分布  q(x) は、サンプルデータ  X^{n} から求めたMAP推定量  \hat{w}_{ML}確率モデル  p(X|w) によって次式のように推定されます。
$$ \begin{eqnarray} q(X) \approx p( X | \hat{w}_{MAP} ) \end{eqnarray} $$



4. ベイズ推定(Bayse)

4.1 パラメータ  w の平均対数尤度  K

ベイズの定理: $$ \begin{eqnarray} p(w | X^{n}) = \frac{ p(X^{n}|w) \varphi(w) }{ \int_{W} p(X^{n}|w) \varphi(w) dw} \end{eqnarray} $$

において、以下で定義される平均対数尤度  K は、 X に対して期待値をとった(  X のばらつきによる偶然誤差を排除した)上での、確率モデル  p(X|w) の尤もらしさを示します。

$$ \begin{eqnarray} K(w) & := & \mathbb{E}_{X} \left[ \log p(X | w) \right] \tag{期待対数尤度} \\ \\ K_n(w) & := & \frac{1}{n} \sum_{i=1}^{n} \log p(X_i | w) \tag{経験対数尤度} \end{eqnarray} $$

パラメータの推定量を求めるため、平均対数尤度を符号反転したものを、損失関数  L として定義します。

$$ \begin{eqnarray} L(w) & := & - K(w) = \mathbb{E}_{X} \left[ - \log p(X | w) \right] \tag{期待損失関数} \\ \\ L_n(w) & := & - K_{n}(w) = \frac{1}{n} \sum_{i=1}^{n} \left\{ - \log p(X_i | w) \tag{経験損失関数} \right\} \end{eqnarray} $$

4.2 確率分布  q(x) の推定(予測分布)

 ベイズ推定では、損失関数が最小となるようなパラメータ  \hat{w}_{0} を選び、これを確率モデル  p(X|w) と事前分布  \varphi(w) によって求められる事後分布  p(w|X^{n}) の平均値の推定量とみなします。ここで、最尤推定やMAP推定では、「確率モデル  p(X|w) に対して与えるパラメータ  w はある1つの値に決まる」と仮定していますが、ベイズ推定では「パラメータ  w は事前分布  \varphi(w) によって確率的に変動する値(=確率変数)である」と仮定していることに注意します。推定量  \hat{w}_{0}ベイズ定量*1と呼びます。

$$ \begin{eqnarray} \hat{w}_{0} := \underset{w \in W}{\rm argmax} ~ L_{n}(w) \end{eqnarray} $$

 推定量  \hat{ w }_{0} は「事後分布において  w がとり得ると期待される(尤もらしい)値」と考えられます。しかし、ベイズ推定において、確率変数  X の真の分布  q(x) を、 p( X | \hat{w}_{MAP} ) によって推定することはしません。その代わり、推定量  \hat{ w }_{0} を平均にもつ事後分布  p(w|X^{n}) から求められる予測分布  p^{*}(X) を用いて真の分布  q(x) を推定します。

$$ \begin{eqnarray} q(X) \approx p^{*}(X) & := & \int_{W} p(X|w)p(w|X^{n}) dw \end{eqnarray} $$


https://www.researchgate.net/profile/Eric_Van_de_Weg/publication/226192090/figure/download/fig1/AS:393798262771713@1470900230184/Probability-distributions-representing-the-principle-of-Bayesian-analysis-The-posterior.png

 上図は、確率モデル・事前分布・事後分布がすべて正規分布に従う場合の各分布を表す。Prior(事前分布)の「ばらつき」に対して、Likelihood(尤度, 最尤推定で用いる)と Posterior(事後分布, MAP推定で用いる)は、より小さい「ばらつき」をもつ。

5. 参考図書

ベイズ統計の理論と方法

ベイズ統計の理論と方法

*1:他の言い方もあるようです。

確率変数に関する4つの「収束概念」

 確率変数の収束(Convergence of random variables) についてまとめてみました。数理統計学において極限現象の解析は特に重要で、漸近理論、極値理論、大偏差理論などで収束概念は欠かせません。
 解析学において、数列の極限から関数の極限を定義したように、統計学(確率論)においても確率変数列の極限から確率変数の極限を定義します。

https://upload.wikimedia.org/wikipedia/commons/4/49/Notions-of-convergence-in-probability-theory.jpg

1. 確率収束 : convergence in probability

確率変数列  U_n := \left\{ U_1, U_2, \dots \right\}に対して、任意の(どんなに小さな)  {\varepsilon}>0について、 $$ \begin{eqnarray} \lim_{n \to \infty} P(~|U_n-U|{\geq}{\varepsilon}~) = 0 \end{eqnarray} $$ が成り立つとき、「確率変数列 \left\{ U_n \right\} は、確率変数  U確率収束する」といい、これを $$ \begin{eqnarray} U_n \underset{~~~~~p}{\to} U \end{eqnarray} $$ と表記する。


 

2. p次平均収束 : convergence in mean

確率変数列  U_n := \left\{ U_1, U_2, \dots \right\}に対して、 $$ \begin{eqnarray} \lim_{n \to \infty} E \left[~{|U_n-U|}^{~p} \right] = 0 \end{eqnarray} $$ が成り立つとき,「確率変数列 \left\{ U_n \right\} は、確率変数  Up次平均収束する」という。  
 

3. 分布収束(法則収束): convergence in distribution (law)

確率変数列  U_n := \left\{ U_1, U_2, \dots \right\}に対して、確率変数  U の分布関数  F_U(x) の任意の連続点  x で、 $$ \begin{eqnarray} \lim_{n \to \infty} P( U_n \leq x ) & = & P( U \leq x) \\ & = & F_U(x) \end{eqnarray} $$ が成り立つとき、「確率変数列 \left\{ U_n \right\} は、確率変数  U分布収束(法則収束)する」といい、これを、 $$ \begin{eqnarray} U_n \underset{~~~~~d}{\to} U \end{eqnarray} $$ と表記する。  
 

4. 概収束 : almost sure convergence

確率変数列  U_n := \left\{ U_1, U_2, \dots \right\}に対して、 $$ \begin{eqnarray} P\left( \omega~{\large |}~{\lim_{n \to \infty} \left| U_n(\omega)-U( \omega ) \right| =0} \right) = 1 \end{eqnarray} $$

が成り立つとき,「確率変数列 \left\{ U_n \right\} は、確率変数  U概収束する」といい,これを、 $$ \begin{eqnarray} U_n \to U ~~~ a.s. \end{eqnarray} $$ と表記する。          

参考

現代数理統計学の基礎 (共立講座 数学の魅力)

現代数理統計学の基礎 (共立講座 数学の魅力)

現代数理統計学 (創文社現代経済学選書)

現代数理統計学 (創文社現代経済学選書)

入門確率過程

入門確率過程

http://kristen-ressurs.no/images/Fysikk/Hilbert-space.jpg http://radhakrishna.typepad.com/.a/6a00d83453b94569e2014e87c56ff9970d-pi