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

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

19年8月の振り返り

8月はずいぶん相場が乱高下してました。まだ当分不安定な状況が続くんですかね。

中断してしまっていた機械学習の勉強をそろそろ再開します。思い出すのが大変なんですよね。思い出すために自分の書いた記事を読み直すんですが、書いたときはわかりやすいと思っていた内容が、読み直すとけっこうひどい。内容が難しくなるにつれて理解の自信も落ちるので、断定的な表現を避けた文章を書いてしまっているか、無意識に曖昧にしてしまっているのかもしれない。いや、まだ読み直していないので想像ですけど。

19年8月の実績
Web収入 投信 前月比評価損益
724,572円 -660,716円 63,856円
19年の累積
Web収入 投信評価損益 個別株
5,141,724円 1,014,752円 705,618円 6,862,094円

19年7月の振り返り

7月はWeb収入も投信も調子が良かった。最近はやりたいことが多くて機械学習に時間が割けていない。おそらく機械学習の記事を期待してくださっているであろう36名の読者の方へ、2回連続でお金の記事を流してしまうことに罪悪感を感じる。

19年7月の実績
Web収入 投信 前月比評価損益
758,917円 367,059円 1,125,976円
19年の累積
Web収入 投信評価損益 個別株
4,417,152円 1,675,468円 705,618円 6,798,238円

19年6月の振り返り

前月はずいぶん相場が下がったけど、6月はちょっと持ち直した。今日7月1日は日経平均も大幅に上昇して7月は幸先の良いスタート。といっても、昨年の-300万が大きすぎて、まだ含み損なんですけど。プラ転まであと一息。 Web収入は大型連休のあった先月よりなぜか上昇。単純計算だと今年のWeb収入は700万いくかどうかってところでしょうか。

19年6月の実績
Web収入 投信 前月比評価損益
622,752円 340,576円 963,328円
19年の累積
Web収入 投信評価損益 個別株
3,658,235円 1,308,409円 705,618円 5,672,262円

ガウス過程による回帰(3)実験結果

概要

今回は、以下の記事を経て導出した結果


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

を使って、ガウス過程による回帰を実験してみました。

www.iwanttobeacat.com

www.iwanttobeacat.com

www.iwanttobeacat.com

ガウス過程による回帰の実験結果

早速結果です。カーネル関数はガウスカーネルを使いました。また、予測分布p(t_{N+1}| \mathbf{t})も求めてみました。 カーネルを用いた回帰でも、カーネル回帰分析(2)実験結果とはずいぶん傾向が違いますね。今回の訓練データの場合では、ガウスカーネルでなくても、多項式モデルから\mathbf{K}を計算してもちゃんと予測分布が求められました。事前分布のパラメータ\alphaの定数倍を微調整する必要がありましたが。ちなみに、線形回帰をベイズ推定で解く(2)予測分布をプロットで求めた予測分布は、線形モデルのパラメータ\mathbf{w}の事前分布を設定し、訓練データから事後分布を求め、\mathbf{w}について積分することで予測分布を求めていました。ガウス過程による回帰も、出発点はパラメータ\mathbf{w}の事前分布を設定することでしたが、考え方がいろいろあって面白いですね。

ガウスカーネルの分散を小さくしてみたら面白い結果になりました。 ガウスカーネルの分散が小さいということは、隣の訓練データの影響が小さくなるのでいびつな回帰曲線になってしまうんですね。

サンプリング

ガウス過程で実験してみたようにサンプリングしてみました。



左がガウス過程の記事でサンプリングした結果。yの期待値は0でしたから、0付近をうろうろした結果となっています。一方で平均値が式(1)で求められた右のサンプリングはどれも訓練データである正弦波を再現するようなサンプリングとなりました。

今回のコードです。

# ガウス過程による回帰
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm


def k(x0, x1):
    sigma = 0.3
    return np.exp(-(x0 - x1) ** 2 / (2 * sigma ** 2))
    # ret = 0
    # for i in range(5):
    #  ret += x0**i * x1**i
    # return ret/0.01


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


# p(t_N+1|t)
def pt(x1, y1):
    k_t = np.vectorize(k)(x.T, np.repeat(x1, N))
    mu = k_t @ a
    v = k(x1, x1) + 1.0/BETA - k_t @ Cn_inv @ k_t
    return np.exp(-(y1 - mu) ** 2 / (2 * v)) / np.sqrt(2*np.pi*v)


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

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

# 係数aを求める
BETA = 5.0
Cn_inv = np.linalg.inv(gram_matrix(x) + 1.0/BETA*np.eye(N))
a = Cn_inv @ t


# 新たな入力x2に対する予測値yを求める
x2 = np.linspace(0, 1, 100)
y = np.empty(100)
for i in range(x2.size):
    y[i] = np.vectorize(k)(x.T, np.repeat(x2[i], N)) @ a


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

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

ガウス過程による回帰(2)

概要

前回の記事の続きです。

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の値がわかれば1つに定まりますので、式(2)の指数部だけに注目してみます。多変量正規分布の確率密度関数は、多変量正規分布の式(5)ですから、p(t_{N+1},\mathbf{t})の指数部は、ガウス過程による回帰(1)の式(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})の指数部はガウス過程による回帰(1)の式(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}

です。

ここで、ガウス過程による回帰(1)の式(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}の条件の表記は省略しています

ガウス過程による回帰(1)

概要

以下の記事で線形回帰において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}

です。ここで誤差\epsilonN(\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の定数倍はカーネル関数の中に含まれているとします。)

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


\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}の条件は省略しています