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

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

19年11月の振り返り

不労所得

11月も相場は好調。昨年は300万くらいマイナスだったので、１年かけて回復したという感じ。もともとはWeb収入の減少分を投資でカバーという作戦だったので、再スタートラインに立ったというところ。

19年11月の実績
Web収入 投信 前月比評価損益 計
630,362円 569,187円 1,199,549円
19年の累積
Web収入 投信評価損益 個別株 計
7,160,792‬円 2,811,288円 705,618円 10,677,698‬円

Woodburyの公式

線形代数・ベクトル

以下の式

(\mathbf{A}+\mathbf{B}\mathbf{C}\mathbf{D})^{-1} = \mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{B}(\mathbf{C}^{-1}+\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1}\mathbf{D}\mathbf{A}^{-1} \tag{1}

は、Woodburyの公式と呼ばれます。左辺の計算より、右辺の計算が楽な場合において用いられるようです。

以下、導出を確認しました。

まず、(\mathbf{I}+\mathbf{P})^{-1}を考えます。

\begin{eqnarray*} (\mathbf{I}+\mathbf{P})^{-1} &=& (\mathbf{I}+\mathbf{P})^{-1}(\mathbf{I}+\mathbf{P}-\mathbf{P}) \tag{3} \\ &=& (\mathbf{I}+\mathbf{P})^{-1} \left\{(\mathbf{I}+\mathbf{P})-\mathbf{P} \right\} \tag{4} \\ &=& \mathbf{I}-(\mathbf{I}+\mathbf{P})^{-1}\mathbf{P} \tag{5} \end{eqnarray*}

です。次に \mathbf{P}+\mathbf{P}\mathbf{Q}\mathbf{P}を考えます。

\begin{eqnarray*} \mathbf{P}+\mathbf{P}\mathbf{Q}\mathbf{P} &=& \mathbf{P}(\mathbf{I}+\mathbf{Q}\mathbf{P}) \tag{6} \\ &=& (\mathbf{I}+\mathbf{P}\mathbf{Q})\mathbf{P} \tag{7} \\ \end{eqnarray*}

ですから、式(6)、(7)より

\mathbf{P} = (\mathbf{I}+\mathbf{P}\mathbf{Q}) \mathbf{P} (\mathbf{I}+\mathbf{P}\mathbf{Q})^{-1} \tag{8}

(\mathbf{I}+\mathbf{P}\mathbf{Q})^{-1} \mathbf{P} = \mathbf{P} (\mathbf{I}+\mathbf{P}\mathbf{Q})^{-1} \tag{9}

です。

ここからが本題で、(\mathbf{A}+\mathbf{B}\mathbf{C}\mathbf{D})^{-1}を変形していきます。

\begin{eqnarray*} (\mathbf{A}+\mathbf{B}\mathbf{C}\mathbf{D})^{-1} &=& \left\{ \mathbf{A}(\mathbf{I}+\mathbf{A}^{-1}\mathbf{B}\mathbf{C}\mathbf{D})\right\}^{-1} \tag{10} \end{eqnarray*}

で、逆行列の定理の式(1)より、

式(10)= (\mathbf{I}+\mathbf{A}^{-1}\mathbf{B}\mathbf{C}\mathbf{D})^{-1}\mathbf{A}^{-1} \tag{11}

です。ここで式(5)を使えば、

式(11)= \left\{ \mathbf{I}-(\mathbf{I}+\mathbf{A}^{-1}\mathbf{B}\mathbf{C}\mathbf{D})^{-1}\mathbf{A}^{-1}\mathbf{B}\mathbf{C}\mathbf{D} \right\}\mathbf{A}^{-1} \tag{12}

となります。これを展開すれば

式(12)= \mathbf{A}^{-1} - ( \mathbf{I}+\mathbf{A}^{-1}\mathbf{B}\mathbf{C}\mathbf{D})^{-1}\mathbf{A}^{-1}\mathbf{B}\mathbf{C}\mathbf{D} \mathbf{A}^{-1} \tag{13}

です。ここで式(9)の関係を2回繰り返し使って、

\begin{eqnarray*} 式(13) &=& \mathbf{A}^{-1} - \mathbf{A}^{-1}( \mathbf{I}+\mathbf{B}\mathbf{C}\mathbf{D}\mathbf{A}^{-1})^{-1}\mathbf{B}\mathbf{C}\mathbf{D} \mathbf{A}^{-1} \tag{14} \\ &=& \mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{B}( \mathbf{I}+\mathbf{C}\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1}\mathbf{C}\mathbf{D} \mathbf{A}^{-1} \tag{15} \end{eqnarray*}

となります。ここで式(15)の( \mathbf{I}+\mathbf{C}\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1}\mathbf{C}について変形していくと、

\begin{eqnarray*} ( \mathbf{I}+\mathbf{C}\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1}\mathbf{C} &=& ( \mathbf{C}\mathbf{C}^{-1}+\mathbf{C}\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1}\mathbf{C} \tag{16} \\ &=& \left\{ \mathbf{C} (\mathbf{C}^{-1}+\mathbf{D}\mathbf{A}^{-1}\mathbf{B}) \right\} ^{-1}\mathbf{C} \tag{17} \\ &=& (\mathbf{C}^{-1}+\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1} \mathbf{C}^{-1}\mathbf{C} \tag{18} \\ &=& (\mathbf{C}^{-1}+\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1} \tag{19} \end{eqnarray*}

となります。途中で逆行列の定理の式(1)を使っています。そして式(19)と式(15)より、

(\mathbf{A}+\mathbf{B}\mathbf{C}\mathbf{D})^{-1} = \mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{B}(\mathbf{C}^{-1}+\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1}\mathbf{D}\mathbf{A}^{-1} \tag{20}

となります。

ガウス過程による分類（３）

カーネル法

概要

前回の以下の記事の続きです。

www.iwanttobeacat.com

前回までのおさらい

訓練データ\mathbf{x}=\{\mathbf{x}_{1},\cdots,\mathbf{x}_{N} \}、\mathbf{t}_{N}=(t_1,\cdots,t_{N})^{T}が与えられたとき、新たな入力\mathbf{x}_{N+1}に対するt_{N+1}=1の確率、すなわちp(t_{N+1}=1 | \mathbf{t}_{N},\mathbf{x},\mathbf{x}_{N+1} )を求めることが目標です。そしてこれは、いくつかの式変形によって以下のように表せました。

\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} )= \int \underline{ p(t_{N+1}=1 |a_{N+1} ) }p(a_{N+1}|\mathbf{t}_{N}) d a_{N+1} \tag{1}

式(1)下線部はシグモイド関数の出力 \sigma(a_{N+1})ですが、p(a_{N+1}|\mathbf{t}_{N})は解析的に求めらません。これは

\displaystyle p(a_{N+1}|\mathbf{t}_{N}) = \begin{align} \int \underset{A} {\underline{p(a_{N+1}|\mathbf{a}_{N}) }} \, \underset{B} {\underline{p(\mathbf{a}_{N}|\mathbf{t}_{N}) }} d\mathbf{a}_{N} \end{align} \tag{2}

と変形することができました。さらに式(2)下線部Aは

p(a_{N+1}|\mathbf{a}_{N}) = N(a_{N+1}| \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{a}_{N}, c- \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}) \tag{3}

と求めることができました。そして式(2)下線部Bのp(\mathbf{a}_{N}|\mathbf{t}_{N})をラプラス近似で正規分布の形にして式(2)を計算していく方針です。

p(\mathbf{a}_{N}|\mathbf{t}_{N})の計算

ベイズの定理*1より、

\displaystyle p(\mathbf{a}_{N}|\mathbf{t}_{N}) = \frac{p(\mathbf{t}_{N} | \mathbf{a}_{N})p(\mathbf{a}_{N})}{p(\mathbf{t}_{N})} \tag{4}

です。正規化項であるp(\mathbf{t}_{N})を無視すれば、

p(\mathbf{a}_{N}|\mathbf{t}_{N}) \approx p(\mathbf{t}_{N} | \mathbf{a}_{N})p(\mathbf{a}_{N}) \tag{5}

です。ラプラス近似を行うためには、関数の対数をとり、そのヘッセ行列などを求めなくてはなりませんので、

\ln p(\mathbf{a}_{N}|\mathbf{t}_{N}) \approx \ln (p(\mathbf{t}_{N} | \mathbf{a}_{N})p(\mathbf{a}_{N})) \tag{6}

を考えます。式(6)を \Psi(\mathbf{a}_{N})とおいて、

\begin{eqnarray*} \Psi(\mathbf{a}_{N}) &=& \ln ( p(\mathbf{t}_{N} | \mathbf{a}_{N}) p(\mathbf{a}_{N})) \tag{7} \\ &=& \ln p(\mathbf{t}_{N} | \mathbf{a}_{N}) + \ln p(\mathbf{a}_{N}) \tag{8} \end{eqnarray*}

です。ここで、まず式(8)の右辺p(\mathbf{t}_{N} | \mathbf{a}_{N})を考えます。p(t_{n}=1|a_{n})=\sigma(a_{n})、p(t_{n}=0|a_{n})=1-\sigma(a_{n})なので、

\displaystyle p(\mathbf{t}_{N} | \mathbf{a}_{N}) =\prod_{n=1}^{N} \sigma(a_{n})^{t_{n}} \left( 1-\sigma(a_{n}) \right)^{1-t_{n}} \tag{9}

と書けます。\sigma(a)はシグモイド関数で、

\displaystyle \sigma(a) = \frac{1}{1+\exp(-a)} \tag{10}

ですから、

\displaystyle \sigma(a)^{t} \left( 1-\sigma(a) \right)^{1-t} = \left( \frac{1}{1+\exp(-a)} \right)^{t} \left( 1- \frac{1}{1-\exp(-a)}\right)^{1-t} \tag{11}

です。式(11)右辺を変形していけば、

\begin{eqnarray*} &=& \left( \frac{1}{1+\exp(-a)} \right)^{t} \left( \frac{\exp(-a)}{1-\exp(-a)}\right)^{1-t}\tag{12} \\ &=& \frac{1}{ ( 1+\exp(-a) )^{t} } \frac{\exp(-a+at)}{ ( 1+\exp(-a) )^{1-t} } \tag{13}\\ &=& \frac{\exp(-a+at)}{ 1+\exp(-a) } \tag{14}\\ &=& \exp(at) \frac{\exp(-a)}{ 1+\exp(-a) } \tag{15}\\ &=& \exp(at) \frac{1}{ \exp(a)+1 } \tag{16}\\ &=& \exp(at) \sigma(-a) \tag{17}\\ \end{eqnarray*}

と計算できます。従って式(9)は、

\displaystyle p(\mathbf{t}_{N} | \mathbf{a}_{N}) =\prod_{n=1}^{N} \exp(a_{n}t_{n}) \sigma(-a_{n}) \tag{18}

です。

そして、式(8)右辺のもう１つの項p(\mathbf{a}_{N})は、ガウス過程という前提条件で、

p(\mathbf{a}_{N}) = N(\mathbf{a}_{N}|\mathbf{0},\mathbf{C}_{N}) \tag{19}

でした。

以上より、

\displaystyle \Psi(\mathbf{a}_{N}) = \ln \prod_{n=1}^{N} \exp(a_{n}t_{n}) \sigma(-a_{n}) + \ln N(\mathbf{a}_{N}|\mathbf{0},\mathbf{C}_{N}) \tag{20} \\

です。いろいろ計算がややこしくて何をやっているか見失いそうになるんですが、今はp(\mathbf{a}_{N}|\mathbf{t}_{N})をラプラス近似で正規分布の形にもっていこうとしている途中です、、、。

続きは次回 www.iwanttobeacat.com

*1:ベイズの定理の式(2)

ガウス過程による分類（２）

カーネル法

概要

前回の以下の記事の続きです。 www.iwanttobeacat.com

前回の記事までのまとめ

訓練データ\mathbf{x}=\{\mathbf{x}_{1},\cdots,\mathbf{x}_{N} \}、\mathbf{t}_{N}=(t_1,\cdots,t_{N})^{T}が与えられたとき、新たな入力\mathbf{x}_{N+1}に対するt_{N+1}=1の確率、すなわちp(t_{N+1}=1 | \mathbf{t}_{N},\mathbf{x},\mathbf{x}_{N+1} )を求めることが目標です。そしてこれは、いくつかの式変形によって以下のように表せました。

\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} )= \int p(t_{N+1}=1 |a_{N+1} )p(a_{N+1}|\mathbf{t}_{N}) d a_{N+1} \tag{1}

p(a_{N+1}|\mathbf{t}_{N})の計算

式(1)において、p(t_{N+1}=1|a_{N+1})=\sigma(a_{N+1})ですから、式(1)右辺の前半は簡単に求まりますが、p(a_{N+1}|\mathbf{t}_{N})は解析的に求められず、これを何とかしていかなくてはなりません。

確率の周辺化の式を用いて、*1

\displaystyle p(a_{N+1}|\mathbf{t}_{N}) = \int p(a_{N+1},\mathbf{a}_{N}|\mathbf{t}_{N}) d\mathbf{a}_{N} \tag{2}

と変形し、さらにベイズの定理により、*2

\displaystyle 式(2)=\frac{1}{p(\mathbf{t}_{N})} \int p(a_{N+1},\mathbf{a}_{N})p(\mathbf{t}_{N}|a_{N+1},\mathbf{a}_{N}) d\mathbf{a}_{N} \tag{3}

です。そして確率の乗法定理により、*3

\displaystyle 式(3)=\frac{1}{p(\mathbf{t}_{N})} \int p(a_{N+1}|\mathbf{a}_{N})p(\mathbf{a}_{N})p(\mathbf{t}_{N}|a_{N+1},\mathbf{a}_{N}) d\mathbf{a}_{N} \tag{4}

です。ここでa_{N+1}と\mathbf{t}_{N}は独立ですから条件から取り除くことができ、

\displaystyle 式(4)=\frac{1}{p(\mathbf{t}_{N})} \int p(a_{N+1}|\mathbf{a}_{N})p(\mathbf{a}_{N})p(\mathbf{t}_{N}| \mathbf{a}_{N}) d\mathbf{a}_{N} \tag{5}

です。そしてベイズの定理の式(1)より、 p(\mathbf{a}_{N}) p(\mathbf{t}_{N} | \mathbf{a}_{N}) = p(\mathbf{t}_{N})p(\mathbf{a}_{N} | \mathbf{t}_{N}) ですから、以上より

\displaystyle \displaystyle p(a_{N+1}|\mathbf{t}_{N}) = \int p(a_{N+1}|\mathbf{a}_{N}) p(\mathbf{a}_{N}|\mathbf{t}_{N}) d\mathbf{a}_{N} \tag{6}

となります。

p(a_{N+1}|\mathbf{a}_{N})の計算

ここでまずは式(6)のp(a_{N+1}|\mathbf{a}_{N})を考えます。

aはガウス過程という前提条件でした。\mathbf{y}がガウス過程に従うとき、 p(\mathbf{y}) = N(\mathbf{y}|\mathbf{0},\alpha^{-1}\mathbf{K}) と書けましたから*4、a_{N+1}もこの形で書くことができるはずです。ここで\alphaの定数倍をカーネル関数の中に含めれば

p(a_{N+1}) = N(a_{N+1}|\mathbf{0},\mathbf{K}) \tag{7}

と書けます。今考えているのは分類問題で離散値データですから、ガウス過程による回帰（１）の式(1)ようにノイズはありませんが、共分散行列の正定値性の保証のため*5、

C(x_{n},x_{m}) = k(x_{n},x_{m}) + \nu \delta_{nm} \tag{8}

のようにノイズの項を加えたものを考えます。これが正定値行列であることが、後々最適解を持つことの保証になります。行列の各要素が式(8)で与えられる共分散行列\mathbf{C}_{N+1}を用いれば、

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

と書けます。すると、ガウス過程による回帰（１）、ガウス過程による回帰（２）でp(t_{N+1}|\mathbf{t}_{N})を計算したのと同様の手順で、

p(a_{N+1}|\mathbf{a}_{N}) = N(a_{N+1}| \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{a}_{N}, c- \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}) \tag{10}

と計算できます。

残りはp(\mathbf{a}_{N}|\mathbf{t}_{N})

次に、式(6)のp(\mathbf{a}_{N}|\mathbf{t}_{N})の計算です。計算の方針はラプラス近似を使って正規分布の形にもっていくことです。そして式(10)とあわせて2つの正規分布の計算とします。

続きは次回 www.iwanttobeacat.com

*1:条件付き確率、同時確率、周辺確率の式(10)

*2:ベイズの定理の式(2)

*3:条件付き確率、同時確率、周辺確率の式(2)

*4:ガウス過程の式(15)

*5:分散は0以上であることに注意し、共分散行列の対角成分が正の場合の２次形式を計算してみれば、常に正になることがわかる。

19年10月の振り返り

不労所得

ようやく投信がプラ転！一番ひどいときは250万くらい含み損がありましたが、なんとか復活してくれました。マイナスだったものが少しプラスになっただけで、決して得をしてるわけじゃないんですが、ずいぶん儲かったような錯覚をしてしまう。とりあえず今年はなんとかプラスをキープして終えて欲しいところ。

19年10月の実績
Web収入 投信 前月比評価損益 計
688,326円 769,582円 1,457,908円
19年の累積
Web収入 投信評価損益 個別株 計
6,530,430円 2,242,101円 705,618円 9,478,149円

ガウス過程による分類（１）

カーネル法

概要

ガウス過程を用いた分類を考えていきます。2クラスの分類なら、確率的生成モデル - シグモイド関数で導出したシグモイド関数を用いて

y=\sigma(a(\mathbf{x})) \tag{1}

とモデル化できます。線形識別ではa(\mathbf{x})の部分は\mathbf{w}^{T}\mathbf{x}のように、係数\mathbf{w}の線形結合を考えていましたが、これを今回はガウス過程a(\mathbf{x})によってモデル化します。

a(\mathbf{x})はガウス過程なので、仮に\mathbf{x}が1次元の入力なら下図左のようなサンプルが得られます。これをシグモイド関数に通したものが下図右で、出力を確率として扱うことができるようになる、というのは確率的生成モデル - シグモイド関数で見たとおりです。縦軸で0.5を閾値にどちらのクラスかを識別します。

2次元入力なら下図のようなイメージです。左がガウス過程から得られた分布、真ん中がそれをシグモイド関数に通したものです。これを識別の結果でよく使う等高線図にしたものが右のグラフです。ガウス過程による分類では、これを訓練データをもとにして期待する識別境界となるようにしていきます。

新たな入力に対する確率分布

さて、目標は、訓練データ\mathbf{x}=\{\mathbf{x}_{1},\cdots,\mathbf{x}_{N} \}、\mathbf{t}_{N}=(t_1,\cdots,t_{N})^{T}が与えられたとき、新たな入力\mathbf{x}_{N+1}に対するt_{N+1}=1の確率、すなわちp(t_{N+1}=1 | \mathbf{t}_{N},\mathbf{x},\mathbf{x}_{N+1} )を求めることです。なお2クラス分類を考えていますので、t\in\{0,1\}とすれば、t_{N+1}=1の確率が求まればt_{N+1}=0の確率も当然求まります。また、以降は入力の条件を省略してp(t_{N+1}=1 | \mathbf{t}_{N} )と表記します。

確率の周辺化により、*1

\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} ) = \int p(t_{N+1}=1 ,a_{N+1} | \mathbf{t}_{N} ) d a_{N+1} \tag{2}

と変形できます。ここでa_{N+1}=a(\mathbf{x}_{N+1})です。つまり新たな入力点に対応したガウス過程の出力値です。2次元入力の上のグラフ左図でいえば、新たな入力座標に対するZ軸の値です。

さらに条件付き同時確率の式変形により式(2)の右辺は、*2

\displaystyle \int p(t_{N+1}=1 ,a_{N+1} | \mathbf{t}_{N} ) d a_{N+1}= \int p(t_{N+1}=1 |a_{N+1} , \mathbf{t}_{N} )p(a_{N+1}|\mathbf{t}_{N}) d a_{N+1} \tag{3}

となります。ここで、p(t_{N+1}=1)は、a_{N+1}が与えられればp(t_{N+1}=1)=\sigma(a_{N+1})で決まりますから、\mathbf{t}_{N}には依存しません。従って、

\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} )= \int p(t_{N+1}=1 |a_{N+1} )p(a_{N+1}|\mathbf{t}_{N}) d a_{N+1} \tag{4}

となります。この辺りの式変形は線形回帰をベイズ推定で解く（１）予測分布の導出とけっこう似ていますね。

で、この式(4)をなんとかして求めていくのですが、けっこう長い道のりになります、、、。全6回にわたって計算が続きます。 ちなみに最終結果は次のようになります。

\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} )\simeq \sigma \left( \frac{\mathbf{k}^{T}\mathbf{C}_{N}(\mathbf{t}_{N}-\boldsymbol {\sigma}_{N})}{1+\displaystyle \frac{\pi}{8}\left\{ c-\mathbf{k}^{T}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{k} \right\} } \right) \tag{5}

次回： www.iwanttobeacat.com

最終的な実験結果はこちらです

www.iwanttobeacat.com

*1:条件付き確率、同時確率、周辺確率の式(10)

*2:条件付き確率、同時確率、周辺確率の式(8)