2019-06-22

ガウス過程による回帰（２）

概要

前回の記事の続きです。

www.iwanttobeacat.com

訓練データ $(\mathbf{x},\mathbf{t})$ が得られているとき、新たな入力 $\mathbf{x}_{N+1}$ に対する $t_{N+1}$ を求めること、つまり $p(t_{N+1}| \mathbf{t})$ を求めることが目標です。*1　先回は、

$\displaystyle p(t_{N+1}| \mathbf{t}) = \frac{p(t_{N+1},\mathbf{t})}{p(\mathbf{t})} \tag{1}$

における $p(t_{N+1},\mathbf{t})$ を計算したところまでです。今回は式(1)を計算していきます。

計算・・・

計算の方針は、 $p(\mathbf{a},\mathbf{b})$ が正規分布であるとき、 $p(\mathbf{a}|\mathbf{b})$ もまた正規分布であるという性質を用いることです。この性質より、 $p(t_{N+1},\mathbf{t})$ は正規分布でしたので、 $p(t_{N+1}| \mathbf{t})$ も正規分布となります。つまり、式(1)は、

$\displaystyle p(t_{N+1}| \mathbf{t}) = \frac{1}{\sqrt{2\pi \sigma^{2}} }\exp \left\{ -\frac{(t_{N+1}-\mu)^{2}}{2\sigma^{2}} \right\} \tag{2}$

の形で書けるはずです。正規分布の形は、 $\mu$ と $\sigma$ の値がわかれば１つに定まりますので、式(2)の指数部だけに注目してみます。多変量正規分布の確率密度関数は、多変量正規分布の式(5)ですから、 $p(t_{N+1},\mathbf{t})$ の指数部は、ガウス過程による回帰（１）の式(9)より、

$\displaystyle \exp \left\{ - \frac{1}{2}(\mathbf{t},t_{N+1})^{T} \mathbf{C}_{N+1}^{-1} (\mathbf{t},t_{N+1}) \right\} \tag{3}$

です。そして $p(\mathbf{t})$ の指数部はガウス過程による回帰（１）の式(8)より、

$\displaystyle \exp \left( - \frac{1}{2}\mathbf{t}^{T} \mathbf{C}_{N}^{-1} \mathbf{t} \right) \tag{4}$

です。したがって、 $p(t_{N+1}| \mathbf{t})$ の指数部は、

$\displaystyle \exp \left\{ - \frac{1}{2}(\mathbf{t},t_{N+1})^{T} \mathbf{C}_{N+1}^{-1} (\mathbf{t},t_{N+1}) +\frac{1}{2}\mathbf{t}^{T} \mathbf{C}_{N}^{-1} \mathbf{t} \right\} \tag{5}$

となります。ここで、

$\mathbf{C}_{N+1}^{-1} = \left( \begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \\ \end{array} \right) \tag{6}$

とおいて、式(5)の $\exp$ の中を展開していくと、

$\displaystyle -\frac{\mathbf{D}}{2} \{ t_{N+1}^{2} + \mathbf{D}^{-1}(\mathbf{t}^{T}\mathbf{B}+\mathbf{C}\mathbf{t})t_{N+1} + const \} \tag{7}$

となります。これを、 $\displaystyle -\frac{(t_{N+1}-\mu)^{2}}{2\sigma^{2}}$ の形にすれば $\mu$ と $\sigma$ が求まります。式(7)を平方完成すれば、

$\displaystyle -\frac{\mathbf{D}}{2} \left\{ t_{N+1} + \frac{1}{2}\mathbf{D}^{-1}( \mathbf{t}^{T}\mathbf{B}+\mathbf{C}\mathbf{t}) \right\}^{2} + const \tag{8}$

となります。よって、

$\displaystyle \mu = -\frac{1}{2} \mathbf{D}^{-1} (\mathbf{t}^{T}\mathbf{B}+\mathbf{C}\mathbf{t} ) \tag{9}$

$\displaystyle \sigma^{2} = \mathbf{D}^{-1} \tag{10}$

です。

ここで、ガウス過程による回帰（１）の式(10)より、

$\mathbf{C}_{N+1} = \left( \begin{array}{cc} \mathbf{C}_{N} & \mathbf{k} \\ \mathbf{k}^{T} & c \\ \end{array} \right) \tag{11}$

でしたので、ブロック行列の逆行列の公式より、

$\mathbf{C}_{N+1}^{-1} = \left( \begin{array}{cc} \mathbf{C}_{N}^{-1} + S^{-1}\mathbf{C}_{N}^{-1}\mathbf{k}\mathbf{k}^{T}\mathbf{C}_{N}^{-1} & -S^{-1}\mathbf{C}_{N}^{-1}\mathbf{k} \\ -S^{-1}\mathbf{k}^{T}\mathbf{C}_{N}^{-1} & S^{-1} \\ \end{array} \right) \tag{12}$

$S = c-\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k} \tag{13}$

となります。 $S$ はスカラです。式(6)と式(12)より、

$\begin{eqnarray*} \displaystyle \mu &=& -\frac{1}{2} \mathbf{D}^{-1} (\mathbf{t}^{T}\mathbf{B}+\mathbf{C}\mathbf{t} ) \tag{14}\\ \displaystyle &=& -\frac{S}{2} \left( -S^{-1}\mathbf{t}^{T} \mathbf{C}_{N}^{-1}\mathbf{k} - S^{-1}\mathbf{k}^{T}\mathbf{C}_{N}^{-1} \mathbf{t} \right) \tag{15}\\ \displaystyle &=& -\frac{S}{2} \left\{ -S^{-1}\mathbf{t}^{T} \mathbf{C}_{N}^{-1}\mathbf{k} - S^{-1} (\mathbf{C}_{N}^{-1}\mathbf{k})^{T} \mathbf{t} \right\} \tag{16} \end{eqnarray*}$

となります。 $\mathbf{C}_{N}$ は分散共分散行列ですから対称行列です。対称行列の逆行列もまた逆行列であることに注意。ここで $\mathbf{C}_{N}^{-1}\mathbf{k}$ はベクトルですから、

$\mathbf{t}^{T}\mathbf{C}_{N}^{-1}\mathbf{k} = (\mathbf{C}_{N}^{-1}\mathbf{k})^{T} \mathbf{t} \tag{17}$

ですので、

$\displaystyle \mu = \mathbf{t}^{T}\mathbf{C}_{N}^{-1}\mathbf{k} = \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{t} \tag{18}$

です。また、

$\begin{eqnarray*} \displaystyle \sigma^{2} &=& \mathbf{D}^{-1} \tag{19} \\ &=& c- \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k} \tag{20} \end{eqnarray*}$

です。

以上より新たな入力に対する分布 $p(t_{N+1}| \mathbf{t})$ が定まりました。予測値は分布の平均値とすればいいので、式(18)で求められることになります。

まとめ

さて、ここで

$\mathbf{a}= \mathbf{C}_{N}^{-1}\mathbf{t} \tag{21}$

とすれば、予測値 $m(\mathbf{x}_{N+1})$ は、

$\displaystyle m(\mathbf{x}_{N+1}) = \sum_{n=1}^{N} a_{n} k(\mathbf{x}_n,\mathbf{x}_{N+1}) \tag{22}$

と書くことができ、カーネル回帰分析の式(16)と同様に、カーネル関数 $\times$ 定数で表すことができます。なんかいろいろやってきましたけど、正則化最小二乗法から導出しても、確率分布を導入して導出しても、係数の事前分布を設定して導出しても、結局たどり着くところは同じなんですね。

以下の記事で、この結果を用いて回帰の実験をしてみました。

www.iwanttobeacat.com

※今回の計算は手持ちの参考書にはなく、多くを自力で行ったので、ところどころ厳密性を欠いていたりするかもしれません。

*1: $\mathbf{x}$ 、 $\mathbf{x}_{N+1}$ の条件の表記は省略しています

2019-06-09

ガウス過程による回帰（１）

カーネル法

概要

以下の記事で線形回帰において $y$ の事前分布はカーネル関数を使って表すことができ、そしてガウス過程となっていることを確認しました。

www.iwanttobeacat.com

今回は、ガウス過程を使って回帰問題を考えます。つまり、訓練データをもとに新たな入力に対する予測値をガウス過程を使って考えていきます。

ガウス過程による回帰

線形回帰では、係数 $\mathbf w$ を最小二乗法や最尤推定で求めたり、係数の事前分布を設定してMAP推定しましたが、ガウス過程のモデルでは係数 $\mathbf w$ が式から消されており存在しません。ガウス過程で見たように、 $\mathbf y$ は $\mathbf x$ のカーネル関数から求められるガウス過程でした。ここで、新たな入力 $\mathbf{x}_{N+1}$ を加えたとしても、 $\mathbf y$ はやはりガウス過程に従うはずです。このことを使って新たな入力に対する予測値を求めていきます。

訓練データ $t$ には、理論値 $y$ に誤差 $\epsilon$ が加わっているとすると

$t_n = y_n + \epsilon_n \tag{1}$

です。ここで誤差 $\epsilon$ は $N(\epsilon|0,\beta^{-1})$ の正規分布に従うとすれば

$p(t_n|y_n) = N(t_n|y_n, \beta^{-1}) \tag{2}$

と書けます。 $n=1,\cdots,N$ のデータをまとめて表せば、各誤差は互いに独立ですから

$p(\mathbf{t}|\mathbf{y}) = N(\mathbf{t}|\mathbf{y}, \beta^{-1}\mathbf{I}) \tag{3}$

となります。なおここで $p(\mathbf{y})$ はガウス過程で確認したように、

$p(\mathbf{y}) = N(\mathbf{y}|\mathbf{0}, \mathbf{K}) \tag{4}$

です。（ $\alpha$ の定数倍はカーネル関数の中に含まれているとします。）

さて、線形回帰をベイズ推定で解く（１）予測分布の導出でも使いましたが、正規分布に対する以下の公式を使います。

$\begin{eqnarray*} p(\mathbf x) &=& N(\mathbf x | \mathbf \mu, \Lambda^{-1}) \tag{5} \\ p(\mathbf y | \mathbf x) &=& N(\mathbf y | \mathbf A \mathbf x + \mathbf b, \mathbf{L}^{-1}) \tag{6} \end{eqnarray*}$

であるとき、

$p(\mathbf y) = N(\mathbf y | \mathbf A \mathbf \mu + \mathbf b , \mathbf{L}^{-1} + \mathbf A \mathbf \Lambda^{-1} \mathbf A^{T}) \tag{7}$

$\mathbf{\Lambda}^{-1}=\mathbf{K}$ 、 $\mathbf{\mu} = \mathbf{0}$ 、 $\mathbf{A}=\mathbf{I}$ 、 $\mathbf{b}=\mathbf{0}$ 、 $\mathbf{L}^{-1}=\beta^{-1}\mathbf{I}$ に対応させれば、式(3)と式(4)より

$p(\mathbf{t}) = N(\mathbf{t} | \mathbf{0}, \mathbf{C}_N) \tag{8}$

となり、これで訓練データの分布がわかりました。ここで $\mathbf{C}_N$ は $\beta^{-1}\mathbf{I}+\mathbf{K}$ です。 $\mathbf{y}$ の分散が $\mathbf{K}$ で、誤差の分散が $\beta^{-1} \mathbf{I}$ でしたから、 $\mathbf{y}$ と $\epsilon$ の和である $\mathbf{t}$ の分散が $\beta^{-1}\mathbf{I}+\mathbf{K}$ となるのは、分散の加法性からも明らかです。

さて、回帰分析の目的は、訓練データが得られたあと、新たな入力 $\mathbf{x}_{N+1}$ に対する $t_{N+1}$ を予測することです。つまり $p(t_{N+1}|\mathbf t)$ を求めることです。これを求めるためにまず $p(t_{N+1},\mathbf{t})$ を計算します。*1

$p(t_{N+1},\mathbf{t})$ は、 $n=1,\cdots,N+1$ のデータをまとめて考え、式(3)以降と同様の計算より、

$p(t_{N+1},\mathbf{t}) = N(t_{N+1},\mathbf{t} | \mathbf{0}, \mathbf{C}_{N+1}) \tag{9}$

となります。ここで

$\mathbf{C}_{N+1} = \left( \begin{array}{cc} \mathbf{C}_{N} & \mathbf{k} \\ \mathbf{k}^{T} & c \\ \end{array} \right) \tag{10}$

$\mathbf{k} = ( k(\mathbf{x_{1}},\mathbf{x_{N+1}}) , \cdots , k(\mathbf{x_{N}},\mathbf{x_{N+1}}))^{T} \tag{11}$

$c = k(\mathbf{x_{N+1}},\mathbf{x_{N+1}}) + \beta^{-1} \tag{12}$

です。ぱっと見、式を追うのが面倒ですが、以下のように、 $\mathbf C_{N+1}$ を分けて考えれば、式(10)～式(12)の通りであることがわかります。

冒頭で書いた「新たな入力 $\mathbf{x}_{N+1}$ を加えたとしても、 $\mathbf y$ はやはりガウス過程に従うはず」ということから $p(t_{N+1},\mathbf{t})$ を求めることができました。

さて、条件付き確率の定義より

$\displaystyle p(t_{N+1}| \mathbf{t}) = \frac{p(t_{N+1},\mathbf{t})}{p(\mathbf{t})} \tag{13}$

ですので、ここから $p(t_{N+1}|\mathbf{t})$ を計算していきます。

予測値の導出計算は、以下の記事に続きます。 www.iwanttobeacat.com

最終的なガウス過程による回帰の実験結果は以下です。

www.iwanttobeacat.com

*1:参考書に倣い、表記の単純化のため $\mathbf{x}$ 、 $\mathbf{x}_{N+1}$ の条件は省略しています

2019-06-05

ブロック行列の逆行列

線形代数・ベクトル

ブロック行列

$\mathbf{P} = \left( \begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \\ \end{array} \right) \tag{1}$

の逆行列は、 $\mathbf{A}$ が正則のとき、

$\left( \begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \\ \end{array} \right)^{-1} = \left( \begin{array}{cc} \mathbf{A}^{-1} + \mathbf{A}^{-1}\mathbf{B}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1} & -\mathbf{A}^{-1}\mathbf{B}\mathbf{S}^{-1} \\ -\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1} & \mathbf{S}^{-1} \\ \end{array} \right) \tag{2}$

$\mathbf{S} = \mathbf{D} -\mathbf{C}\mathbf{A}^{-1}\mathbf{B} \tag{3}$

です。

本当にそうなっているか、式(2)の左から $\mathbf{P}$ をかけて確かめてみます。

$\begin{eqnarray*} 1行1列目 &=& \mathbf{A}\mathbf{A}^{-1} + \mathbf{A}\mathbf{A}^{-1}\mathbf{B}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1}-\mathbf{B}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1} \tag{4}\\ &=& \mathbf{I} +\mathbf{B}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1}-\mathbf{B}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1} \tag{5} \\ &=&\mathbf{I} \tag{6} \end{eqnarray*}$

$\begin{eqnarray*} 1行2列目 &=& -\mathbf{A}\mathbf{A}^{-1}\mathbf{B}\mathbf{S}^{-1}+\mathbf{B}\mathbf{S}^{-1} \tag{7} \\ &=& \mathbf{0} \tag{8} \end{eqnarray*}$

$\begin{eqnarray*} 2行1列目 &=& \mathbf{C}\mathbf{A}^{-1} +\mathbf{C}\mathbf{A}^{-1}\mathbf{B}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1}-\mathbf{D}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1} \tag{9} \\ &=& \mathbf{C}\mathbf{A}^{-1} - (\mathbf{D} - \mathbf{C}\mathbf{A}^{-1}\mathbf{B})\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1} \tag{10} \\ &=& \mathbf{C}\mathbf{A}^{-1}-\mathbf{S}\mathbf{S}^{-1}\mathbf{C}\mathbf{A}^{-1} \tag{11} \\ &=& \mathbf{0} \tag{12} \end{eqnarray*}$

$\begin{eqnarray*} 2行2列目 &=& -\mathbf{C}\mathbf{A}^{-1}\mathbf{B}\mathbf{S}^{-1} + \mathbf{D}\mathbf{S}^{-1} \tag{13} \\ &=& (\mathbf{D}-\mathbf{C}\mathbf{A}^{-1}\mathbf{B})\mathbf{S}^{-1} \tag{14} \\ &=& \mathbf{S}\mathbf{S}^{-1} \tag{15} \\ &=& \mathbf{I} \tag{16} \end{eqnarray*}$

ちゃんとなっていました。

$\mathbf{D}$ が正則のときは

$\left( \begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \\ \end{array} \right)^{-1} = \left( \begin{array}{cc} \mathbf{M} & -\mathbf{M}\mathbf{B}\mathbf{D}^{-1} \\ -\mathbf{D}^{-1}\mathbf{C}\mathbf{M} & \mathbf{D}^{-1} + \mathbf{D}^{-1}\mathbf{C}\mathbf{M}\mathbf{B}\mathbf{D}^{-1} \\ \end{array} \right) \tag{17}$

$\mathbf{M} = (\mathbf{A}-\mathbf{B}\mathbf{D}^{-1}\mathbf{C})^{-1} \tag{18}$

となります。

2019-06-04

ガウス過程

カーネル法

概要

今回はガウス過程の導出と、線形回帰における出力 $\mathbf{y}$ の事前分布はカーネル関数を使って表すことができるということの確認です。計算自体は難しくないのですが、ガウス過程という名称から確率過程のように時間軸を意識してしまうと戸惑ってしまいます。

線形回帰における出力 $\mathbf{y}$ の事前分布

まずは線形回帰で何度もやってきた通り、回帰式のモデルを

$y = \mathbf{w}^{T}\mathbf{x} \tag{1}$

と設定します。線形回帰をMAP推定で解くで計算したのと同じく、係数 $\mathbf{w}$ の事前分布を

$p(\mathbf{w}) = N(\mathbf{w}|\mathbf{0},\alpha^{-1}\mathbf{I}) \tag{2}$

とします。平均が0で、各パラメータには相関がなく、全て等しい分散 $\alpha^{-1}$ であるということです。線形回帰をMAP推定で解くの式(5)で $\mathbf m_0=0$ 、 $\mathbf S_0 = \alpha^{-1}\mathbf I$ とすれば同じことです。

ここで最小二乗法の解の導出の式(4)と同じように計画行列 $\mathbf{X}$ を用いて

$\mathbf{y} = \mathbf{X}\mathbf{w} \tag{3}$

とします。 $\mathbf{y}$ は訓練データ $x_1, \cdots , x_N$ に対する値 $y_1, \cdots , y_N$ をまとめたベクトルになっています。冒頭で書いた通り、今回はこの $\mathbf{y}$ の分布を求めることが目標です。

さて、正規分布に従う独立な確率変数の和は、やはり正規分布に従うという性質があります。正規分布に従うと仮定した係数 $\mathbf{w}$ と入力 $\mathbf{x}$ の線形和が $y$ ですから、この性質より $y$ は正規分布に従うことがわかります。*1

今回の目的は $\mathbf{y}$ の分布を求めることでしたから、正規分布に従うということは平均と共分散がわかれば良いということになります。まずは平均を求めますと、

$\begin{eqnarray*} E[\mathbf{y}] &=& E[ \mathbf{X}\mathbf{w}] \tag{4} \\ &=& \mathbf{X} E[\mathbf{w}] \tag{5} \\ &=& \mathbf{0} \tag{6} \end{eqnarray*}$

です。 $\mathbf{w}$ は式(2)で平均 $\mathbf{0}$ としていますから式(5)は $\mathbf{0}$ と計算できます。

次に $\mathbf{y}$ の共分散を求めます。多変量正規分布の式(2)の通り、

$\mathrm{Cov}(y_a,y_b) = E[y_a y_b] - E[y_a]E[y_b] \tag{7}$

です。ここで式(6)より $E [ \mathbf{y} ] = \mathbf{0}$ でしたから、

$\mathrm{Cov}(y_a,y_b) = E[y_a y_b] \tag{8}$

です。共分散行列は多変量正規分布の式(4)ですので、

$\begin{eqnarray*} \mathrm{Cov} [ \mathbf{y} ] &=& E[ \mathbf{y}\mathbf{y}^{T} ] \tag{9} \\ &=& E[ \mathbf{X}\mathbf{w} \mathbf{w}^{T}\mathbf{X}^{T} ] \tag{10} \\ &=&\mathbf{X} E[ \mathbf{w} \mathbf{w}^{T} ] \mathbf{X}^{T} \tag{11} \\ \displaystyle &=& \mathbf{X} \frac{1}{\alpha} \mathbf{I} \mathbf{X}^{T} \tag{12} \\ \displaystyle &=& \frac{1}{\alpha} \mathbf{X} \mathbf{X}^{T} \tag{13} \\ \end{eqnarray*}$

と計算できます。式(9)から式(10)は転置行列の定理を使っています。また、 $E[ \mathbf{w} \mathbf{w}^{T}$ ]は前提条件より $\alpha^{-1}\mathbf{I}$ です。ここでカーネル法による正則化最小二乗法（１）の式(10)より、

$\mathbf{K}= \mathbf{X}\mathbf{X}^{T} \tag{14}$

でした。以上をまとめると、求めたい分布は

$p(\mathbf{y}) = N(\mathbf{y}|\mathbf{0},\alpha^{-1}\mathbf{K}) \tag{15}$

となることがわかります。カーネル関数さえ定まればいきなり $y$ の事前分布が求まってしまうっていうことなんですが、なんか驚きです。係数 $\mathbf w$ の事前分布を設定し、計算の過程で $\mathbf w$ の前提条件から $\mathbf w$ を式からなくしてしまい、 $y$ の事前分布になるという。

ガウス過程

さて、ここでようやくタイトルのガウス過程の話です。ガウス過程とは、任意の個数の入力を選んだとき、それに対する出力が常に多次元正規分布に従うことをいいます。ちょっと説明が正確でないかもしれません、、、Wikiを参考にしました：ガウス過程 - Wikipedia

$p(\mathbf{y})$ は式(15)の通り多次元正規分布ですから、これを満たしています。つまりガウス過程となっています。確率過程のように時間軸が関係しそうな名前ですが、とにかく多数個入力の出力が多次元正規分布に従っていればガウス過程と呼べます。

式(15)を使って、 $\mathbf y$ を出力してみた実験結果です。いくつかサンプリングしてみました。

上のグラフは、5次の多項式モデルを考え、 $\mathbf{K}$ を計算しています。下のグラフはガウスカーネルで $\mathbf{K}$ を計算しました。ガウスという名前がたくさん出てきてややこしいですが、ガウスカーネルを使っていない上の図もガウス過程です。

これらの結果はなめらかな曲線になっていますが、これは $\mathbf{K}$ が実際どのような値になっているかを見れば明らかです。例えば-1～1の範囲で均等に5点 $x$ を取り、この $x$ に対してガウスカーネルで $\mathbf{K}$ を計算すれば、以下のようになります。（ $\sigma=0.3$ としました。）

$\mathbf{K} = \left( \begin{array}{ccccc} 1& 0.249& 0.00387& 3.73E-06& 2.23E-10 \\ 0.249& 1& 0.249& 0.00387& 3.73E-06 \\ 0.00387& 0.249& 1& 0.249& 0.00387 \\ 3.73E-06& 0.00387& 0.249& 1& 0.249 \\ 2.23E-10& 3.73E-06& 0.00387& 0.249& 1 \end{array} \right) \tag{16}$

$x$ の値が近いほうが共分散が高い、つまり同じような値を取る確率が高いことがわかります。また、逆に $x$ の値が遠ければ共分散は0に近く無相関に近くなります。この共分散行列によって $\mathbf{y}$ が生成されるのですから、なめらかな曲線が得られることがわかります。

入力は多次元になっても同じようにサンプリングすることができます。こちらもガウスカーネルを用いました。

ガウス過程を回帰や分類で使う

このガウス過程を回帰や分類問題に使っていきます。

ガウス過程による回帰の記事１つ目：

www.iwanttobeacat.com

回帰の最終的な実験結果はこちらです。

www.iwanttobeacat.com

ガウス過程による分類の記事１つ目

www.iwanttobeacat.com

分類の最終的な実験結果はこちらです。

www.iwanttobeacat.com

今回のコードです。

# ガウス過程からのサンプリング
import matplotlib.pyplot as plt
import numpy as np


# カーネル
def kernel(x1, x2):
    # 5次の多項式モデルを考えたときのカーネル関数（内積）
    ret = 0
    for i in range(5):
        ret += x1**i * x2**i
    return ret
    # ガウスカーネル
    # sigma = 0.3
    # return np.exp(-(x1 - x2) ** 2 / (2 * sigma ** 2))



# グラム行列
def gram_matrix(x):
    gram = np.empty([N, N])
    for i in range(N):
        for j in range(N):
            gram[i, j] = kernel(x[i], x[j])
    return gram


# x
N = 100
x = np.linspace(-1, 1, N)

# α
ALPHA = 1.0

# グラム行列K
K = gram_matrix(x) / ALPHA

for i in range(5):
    y = np.random.multivariate_normal(np.zeros(N), K)
    plt.plot(x, y)


# 結果の表示
plt.xlim(-1, 1)
plt.show()

*1:参考：再生性 - Wikipedia

2019-06-01

19年5月の振り返り

不労所得

5月は連休明けから相場が下がりっぱなし。そのうち落ち着くと放っておいたら月末までずっと下がりっぱなし。いや、どっちみち放っておくことには変わりないんですが。前月比-80万ですけど、体感的にはもっとマイナスかと思ってた。Web収入のほうも10連休があったのに約60万という厳しい結果。去年の春連休は90万くらいだったのに→18年5月の振り返り。どちらもこの辺で下げ止まってくれるといいのですが。Web収入はどうせ半分近く税金で持っていかれるんで、できれば投資のほうで挽回して欲しいのが本音ですけど、どうなるでしょうね。

19年5月の実績

Web収入	投信前月比評価損益	計
578,761円	-805,528円	-226,767円

19年の累積

Web収入	投信評価損益	個別株	計
3,035,483円	967,833円	705,618円	4,708,934円

2019-05-26

カーネル回帰分析（２）実験結果

カーネル法

先回カーネル回帰分析（１）で導出した結果を使って、回帰問題を解いてみたいと思います。せっかくなので回帰曲線だけではなく、予測分布 $p(t|\mathbf{x})$ も求めてみます。

まずは予測分布 $p(t|\mathbf{x})$ を求めるための準備。 $p(t|\mathbf{x})$ は条件付き確率、同時確率の式(1)より、

$\displaystyle p(t|\mathbf{x}) = \frac{p(\mathbf{x},t)}{p(\mathbf{x})} \tag{1}$

です。 $p(\mathbf{x},t)$ は、カーネル回帰分析（１）の式(1)より

$\displaystyle p(\mathbf{x},t) = \frac{1}{N}\sum_{n=1}^{N} f(\mathbf{x}-\mathbf{x}_n,t-t_{n}) \tag{2}$

です。 $p(\mathbf{x})$ は、条件付き確率、同時確率の式(10)より、

$\displaystyle p(\mathbf{x}) = \int p(\mathbf{x},t) dt \tag{3}$

ですから、

$\displaystyle p(\mathbf{x}) = \frac{1}{N}\sum_{n=1}^{N}\int f(\mathbf{x}-\mathbf{x}_n,t-t_{n})dt \tag{3}$

です。ここでカーネル回帰分析（１）の式(13)より、

$\displaystyle g(\mathbf{x}) = \int f(\mathbf{x}, t) dt \tag{4}$

でしたから、

$\displaystyle g(\mathbf{\mathbf{x}-\mathbf{x}_n}) = \int f(\mathbf{x}-\mathbf{x}_n, t) dt \tag{5}$

です。 $t$ は無限区間で積分しますから、 $-t_n$ の項は影響しませんので、

$\displaystyle p(\mathbf{x}) = \frac{1}{N}\sum_{n=1}^{N} g(\mathbf{\mathbf{x}-\mathbf{x}_n}) \tag{6}$

です。これで予測分布 $p(t|\mathbf{x})$ が求められるようになりました。

さて、肝心のカーネル関数はカーネル密度推定法と同じくガウスカーネルを使い、

$\displaystyle g(\mathbf{x}-\mathbf{x}_n) = \frac{1}{(2 \pi h^{2})^{\frac{D}{2}}} \exp \left( -\frac{ \| \mathbf{x}-\mathbf{x}_n \|^{2}}{2h^{2}} \right) \tag{7}$

とします。

ちょっと前置きが長かったですが、実験結果は以下のようになりました。

f:id:opabinia2:20190526230326p:plain まずは回帰曲線。ちゃんと近似できています。次に予測分布。

f:id:opabinia2:20190526230359p:plain

線形回帰をベイズ推定で解く（２）予測分布をプロットでは回帰曲線のモデルを多項式などと設定していましたが、ノンパラメトリックなカーネル回帰分析ではそういう設定がありませんので、同じ予測分布でもパラメトリック手法とは傾向が違って面白いですね。今回ガウスカーネルを使っていますので、予測分布の結果は、訓練データを中心とした正規分布が配置されていると思えばまあこうなるよねって感じですね。

上記の結果は、平滑化パラメータ $h$ を $0.1$ としていますが、 $0.05$ 、 $0.2$ としてみると、それぞれ次のような結果となりました。

f:id:opabinia2:20190526230843p:plain f:id:opabinia2:20190526230832p:plain

このあたりの結果はカーネル密度推定法の実験と同じですね。平滑化パラメータが小さければ訓練データに強く影響される結果となり、大きければ影響が弱くなだらかな分布となります。

今回のコードです。

# カーネル回帰　Nadaraya-Watsonモデル
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm


# g(x-xn)（ガウスカーネル）
def g(x0):
    return np.exp(-(x0-x)**2/(2*(h**2))) / ((2*np.pi*h**2)**(D/2))


# カーネル k(x,x')
def kernel(x0):
    return g(x0)/np.sum(g(x0))


# p(x,t)
def pxt(x0, t0):
    # (x0-x)^2 + (t0-t)^2　は、 ||(x0,t0) - (x,t)||^2 を展開した形
    return np.average(np.exp(-((x0-x)**2 + (t0-t)**2)/(2*(h**2))) / ((2*np.pi*h**2)**(D/2)) )


# p(x)
def px(x0):
    return np.average(g(x0))


# ランダムシードを固定
np.random.seed(0)
# 訓練データ数
N = 10

# ガウスカーネルのパラメータ
h = 0.1
D = 1

# 訓練データ
x = np.linspace(0, 1, N)
t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N)

# 新たな入力x2に対する予測値yを求める
x2 = np.linspace(0, 1, 100)
y = np.empty(100)

# 予測値は訓練データとカーネル関数の線形和
for i in range(100):
    y[i] = t @ kernel(x2[i])

# 結果の表示
plt.xlim(0.0, 1.0)
plt.ylim(-1.5, 1.5)
plt.scatter(x, t)
plt.plot(x2, y)
plt.show()

# 予測分布p(t|x)
X, Y = np.meshgrid(np.linspace(-0, 1, 100), np.linspace(-1.5, 1.5, 100))
Z = np.vectorize(pxt)(X, Y)/np.vectorize(px)(X)
plt.contourf(X, Y, Z, cmap=cm.coolwarm)
plt.scatter(x, t)
plt.xlim(0,1)
plt.ylim(-1.5,1.5)
plt.show()

2019-05-19

カーネル回帰分析（１）

カーネル法

カーネル法による正則化最小二乗法（１）で、正則化最小二乗法の誤差関数からカーネル関数を導きました。今回は、別のアプローチで回帰分析を行ってもやはりカーネル関数が導かれることを見ていきます。

与えられた訓練データ $(\mathbf{x}_n,t_n)$ をもとに、分布 $p(\mathbf{x},t)$ を推定することを考えます。カーネル密度推定法を用い、密度関数（カーネル密度推定法におけるカーネル関数）を $f(\mathbf{x},t)$ とすれば、カーネル密度推定法の式(13)より、

$\displaystyle p(\mathbf{x},t) = \frac{1}{N}\sum_{n=1}^{N} f(\mathbf{x}-\mathbf{x}_n,t-t_{n}) \tag{1}$

と書けます。

さて、求めたいものは新たな入力に対する予測値 $y$ です。ここで $\mathbf{x}$ を固定したときの $t$ の分布 $p(t|\mathbf{x})$ を考えれば、 $p(t|\mathbf{x})$ の平均値を予測値とすれば良さそうです。よって $y$ の値は、

$\begin{eqnarray*} y(\mathbf{x}) &=& \mathbb{E}[t | \mathbf{x}] \tag{2} \\ &=& \int_{-\infty}^{\infty} t p(t|\mathbf{x})dt \tag{3} \end{eqnarray*}$

です。式(3)は確率密度関数に対する期待値の計算を定義どおりに行っているだけです。

条件付き確率、同時確率、周辺確率の式(1)より

$\displaystyle P(X|Y) = \frac{P(X,Y)}{P(Y)} \tag{4}$

ですから、式(3)は

$\displaystyle y(\mathbf{x}) = \frac{1}{p(\mathbf{x})} \int_{-\infty}^{\infty} t p(\mathbf{x},t)dt \tag{5}$

です。ここで条件付き確率、同時確率、周辺確率の式(10)より

$\displaystyle p(\mathbf{x}) = \int_{-\infty}^{\infty}p(\mathbf{x},t)dt \tag{6}$

ですから、式(5)に代入すれば

$\displaystyle y(\mathbf{x}) = \frac{\displaystyle \int_{-\infty}^{\infty} t p(\mathbf{x},t)dt }{\displaystyle \int_{-\infty}^{\infty}p(\mathbf{x},t)dt} \tag{7}$

となります。さらに式(1)を式(7)に代入すれば、

$\displaystyle y(\mathbf{x}) = \frac{\displaystyle \sum_{n=1}^{N} \int_{-\infty}^{\infty} t f(\mathbf{x}-\mathbf{x}_n,t-t_{n})dt }{\displaystyle \sum_{n=1}^{N} \int_{-\infty}^{\infty} f(\mathbf{x}-\mathbf{x}_n,t-t_{n})dt} \tag{8}$

となります。積分の演算の線形性より、 $\sum$ は積分の外に出せます。

さて、ここで任意の $\mathbf{x}$ に対して

$\displaystyle \int t f(\mathbf{x}, t) dt = 0 \tag{9}$

が成り立つとします。式(9)は平均値の計算です。用いた密度関数 $f(\mathbf{x},t)$ を中心0の正規分布と考えれば、とりあえずイメージが湧きそうです。で、式(9)が成り立つとき、 $t$ に $t-t_n$ 、 $\mathbf{x}$ に $\mathbf{x}-\mathbf{x}_n$ を代入すれば、

$\begin{eqnarray*} \displaystyle \int_{-\infty}^{\infty} (t-t_n)f(\mathbf{x}-\mathbf{x}_n, t-t_n)dt &=& 0 \tag{10} \\ \displaystyle \int_{-\infty}^{\infty} tf(\mathbf{x}-\mathbf{x}_n, t-t_n)dt &=& t_n \int_{-\infty}^{\infty} f(\mathbf{x}-\mathbf{x}_n, t-t_n)dt \tag{11} \end{eqnarray*}$

となり、式(11)を式(8)に代入すると、

$\displaystyle y(\mathbf{x}) = \frac{\displaystyle \sum_{n=1}^{N} t_n \int_{-\infty}^{\infty} f(\mathbf{x}-\mathbf{x}_n, t-t_n)dt }{\displaystyle \sum_{n=1}^{N} \int_{-\infty}^{\infty} f(\mathbf{x}-\mathbf{x}_n,t-t_{n})dt} \tag{12}$

となります。ここで

$\displaystyle g(\mathbf{x}) = \int f(\mathbf{x}, t) dt \tag{13}$

と書くことにすれば、

$\displaystyle y(\mathbf{x}) = \frac{\displaystyle \sum_{n=1}^{N} g(\mathbf{x}-\mathbf{x}_n)t_n }{\displaystyle \sum_{n=1}^{N} g(\mathbf{x}-\mathbf{x}_n)} \tag{14}$

と書けます。さらに、カーネル関数を

$k(\mathbf{x},\mathbf{x}^{\prime}) = \frac{\displaystyle g(\mathbf{x}-\mathbf{x}^{\prime}) }{\displaystyle \sum_{n=1}^{N} g(\mathbf{x}-\mathbf{x}_n)} \tag{15}$

とおけば、

$\displaystyle y(\mathbf{x}) = \sum_{n=1}^{N}k(\mathbf{x}, \mathbf{x}_n)t_n \tag{16}$

となり、推定値 $y$ はカーネル関数 $\times$ 定数で表されることがわかります。別のアプローチでも正則化最小二乗法からの導出と同じ結果となりました。で、式(16)を用いた回帰分析をカーネル回帰、またはNadaraya-Watsonモデルと呼ぶようです。

最後のほうはなんか騙されたような気になっちゃったんですが。例えばカーネル関数って内積として定義したんじゃなかったっけ、式(15)はカーネル関数とみなせるのか？とか（参考：カーネル法）。自分なりには、例えば $g(\mathbf{x}-\mathbf{x}^{\prime})$ がカーネル関数なら何らかの内積として表せるはずで、それが $(a_1, \cdots, a_M)^{T}(b_1, \cdots, b_M)$ と書けるなら、式(15)は分母を $h(\mathbf{x})$ と書けば $(\frac{a_1}{\sqrt{h(\mathbf{x})}}, \cdots, \frac{a_M}{\sqrt{h(\mathbf{x})}})^{T}(\frac{b_1}{\sqrt{h(\mathbf{x})}}, \cdots, \frac{b_M}{\sqrt{h(\mathbf{x})}})$ となるはずだから、やっぱりカーネル関数、という理解をしました。合ってるのかな？あんまり自信がない。

まあ式変形はいろいろあるんですが、カーネル密度推定法で分布 $p(\mathbf{x},t)$ を推定し、 $\mathbf{x}$ を固定したときの分布の平均値を予測値とするというだけです。

次回プログラムを書いて実験してみたいと思います。　→　カーネル回帰分析（２）実験結果

この実験結果を見てからのほうが、この記事でやっていることのイメージが湧きやすいかもしれません。

機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100％趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

ガウス過程による回帰（２）

概要

計算・・・

まとめ

ガウス過程による回帰（１）

概要

ガウス過程による回帰

ブロック行列の逆行列

ガウス過程

概要

線形回帰における出力 $\mathbf{y}$ の事前分布

ガウス過程

ガウス過程を回帰や分類で使う

19年5月の振り返り

19年5月の実績

19年の累積

カーネル回帰分析（２）実験結果

カーネル回帰分析（１）

概要

計算・・・

まとめ

概要

ガウス過程による回帰

概要

線形回帰における出力の事前分布

ガウス過程

ガウス過程を回帰や分類で使う

19年5月の実績

19年の累積

線形回帰における出力 $\mathbf{y}$ の事前分布