スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

GWに学んだことシリーズその1 サポートベクトル回帰

サポートベクトル回帰について学んだことを自分用にまとめる. $\bm x$を入力(説明変数),$\bm\phi(\bm x)$を特徴ベクトルとしたとき, 回帰に用いるモデルを \begin{align} \hat y(\bm x) = \bm w^\top \bm \phi(\bm x) + b \end{align} とする.また,N組の説明変数と目的変数 $\bm x_n, $,$y_n$($n = 1,\dots N$)が得られているとする.このとき,評価関数 \begin{align} C\sum^N_{n=1}E_\varepsilon\left(\hat y(\bm x_n)-y_n\right ) + \frac{1}{2}\|\bm w\|^2 \end{align} を最小化するパラメータ$\bm w$,$b$を求めるのがサポートベクトル回帰である.ただし, \begin{align} E_\varepsilon(e) = \begin{cases} 0 & (|e| < \varepsilon) \\ |e| - \varepsilon & (\mathrm{otherwise}) \end{cases} \end{align} である.この最適化問題をいろいろ頑張って書き換えていくと,カーネル関数 \begin{align} k(\bm x, \bm x') = \bm \phi^\top(\bm x)\bm \phi(\bm x') \end{align}を用いてモデルを \begin{align} \hat y(\bm x) = -\sum^N_{n=1}\left(a_n^+ - a_n^-\right)k(\bm x,\bm x_n) + b \end{align}のように書き換えることができて,このパラメータ推定問題は, \begin{align} 0 \leq a_n^+ \leq C\\ 0 \leq a_n^- \leq C \end{align}のもとで \begin{align} L &= \bm a^\top \bm P \bm a + \bm a^\top \bm q\\ \end{align} を最小化する二次計画問題に落ちる. ただし, \begin{align} \bm a &= \begin{bmatrix} a^+_1 & \dots & a^+_N & a_1^- & \dots & a_N^- \end{bmatrix}^\top\\ \bm P &= \begin{bmatrix}\bm K & -\bm K\\ -\bm K& \bm K\end{bmatrix}\\ \bm q &= \begin{bmatrix} y_1 + \varepsilon & \dots & y_N + \varepsilon & -y_1 + \varepsilon & \dots & -y_N + \varepsilon \end{bmatrix}^\top \end{align} である.ここで,$\bm K$はグラム行列と呼ばれて,$(n,m)$成分が$k(\bm x_n, \bm x_m)$で与えられる$N \times N$行列である.
この最適化問題を解くと多くの$n$について$a_n^+ = a_n^- = 0$となり,スパースなモデルが得られる.
このとき, $a_n^+ \neq 0$または$a_n^- \neq 0$が成り立つサンプル$\bm x_n$をサポートベクトルと呼ぶ.

最後に,$b$は$a_n^+ \neq 0$または$a_n^- \neq 0$が成り立つ$n$について \begin{align} \begin{cases} y_n + \varepsilon + \sum^N_{m=1}\left(a_m^+ - a_m^-\right)k(\bm x_n,\bm x_m) & (a_n^+\neq 0)\\ y_n - \varepsilon + \sum^N_{m=1}\left(a_m^+ - a_m^-\right)k(\bm x_n,\bm x_m) & (a_n^-\neq 0) \end{cases} \end{align}を計算して,平均すれば求められる.


以上をpython(numpy, cvxopt, matplotlib)を用いて実装したのが以下のコード.
import numpy as np
import matplotlib.pylab as plt
import cvxopt
import cvxopt.solvers
class SVR:
    sv = None
    coeff = None
    b = None
    def __init__(self, kernel, e, C):
        self.kernel = kernel
        self.epsilon = e
        self.C = C
    def train(self, x, y):
        N = x.shape[0]
        K = self._makeKmat(x)
        P0 = np.vstack([np.hstack([K, -K]), np.hstack([-K, K])])
        q0 = np.hstack([y, -y]) + self.epsilon
        G0 = np.vstack([np.eye(2 * N), -np.eye(2 * N)])
        h0 = np.hstack([np.ones(2 * N) * self.C, np.zeros(2 * N)])
        A0 = np.hstack([np.ones(N), -np.ones(N)])[np.newaxis, :]
        b0 = 0.
        P = cvxopt.matrix(P0)
        A = cvxopt.matrix(A0)
        q = cvxopt.matrix(q0)
        G = cvxopt.matrix(G0)
        h = cvxopt.matrix(h0)
        b = cvxopt.matrix(b0)
        sol = cvxopt.solvers.qp(P, q, G, h, A, b)
        a = np.array(sol['x']).flatten()
        a[a < 1e-5] = 0
        a_tilde = -a[0:N] + a[N:2*N]
        self.sv = x[a_tilde != 0, :]
        self.coeff = a_tilde[a_tilde != 0]
        self.b = 0
        yhat = np.dot(K, a_tilde)
        bhat = y[a[N:2*N]!=0] - self.epsilon - yhat[a[N:2*N]!=0]
        bhat2 = y[a[0:N]!=0] + self.epsilon - yhat[a[0:N]!=0]
        self.b = np.mean(np.hstack([bhat, bhat2]))
        print "number of support vector is", self.coeff.shape[0]
    def _makeKmat(self, x):
        N = x.shape[0]
        K = np.zeros((N, N))
        for i in range(N):
            K[:,i] = self.kernel.get_value(x, x[i, :])
        return K
    def get_value(self, x):
        N = x.shape[0]
        k = np.zeros((N, self.coeff.shape[0]))
        for i in range(N):
            k[i,:] = self.kernel.get_value(x[i, :], self.sv)
        yhat = np.dot(k, self.coeff) + self.b
        return yhat
    
class gaussianKernel:
    def __init__(self, sigma = None):
        if not sigma:
            self.sigma = 1
        else:
            self. sigma = sigma
    def get_value(self, x1, x2):
       x =  x1 - x2
       if x.ndim == 1:
           e2 = np. inner(x, x)
       else:
           e2 = np.sum(x * x, 1)
       val = np.exp(-e2/(2 * self.sigma**2))
       return val

if __name__ == '__main__':
    kernel = gaussianKernel()
    svr = SVR(kernel, 0.1, 100)
    N = 200
    X = np.random.normal(size = (N, 2))
    y = np.sin(X[:,0]) + np.sin(2 * X[:, 1]) + 5 
    svr.train(X, y + 1e-3 * np.random.normal(size = N))
    yhat = svr.get_value(X)
    plt.plot(y)
    plt.plot(yhat)
    plt.show()
スポンサーサイト

数式

数式をブログに書いてみるテスト.
http://genkuroki.web.fc2.com/
を参考にした.


簡単な使い方:次を HTML ファイルの <head> と </head> のあいだに挿入する。それだけで LaTeX 方式で数式を書けるようになる。

<script type="text/x-mathjax-config">
MathJax.Hub.Config({ tex2jax: { inlineMath: [['$','$'], ["\\(","\\)"]] } });
</script>
<script type="text/javascript"
src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML">
</script>
<meta http-equiv="X-UA-Compatible" CONTENT="IE=EmulateIE7" />



太字だけ追加して,以下のように設定した.

<script type="text/x-mathjax-config">
MathJax.Hub.Config({ tex2jax: { inlineMath: [['$','$'], ["\\(","\\)"]] },
TeX:{Macros:{
bm: ["\\boldsymbol{#1}", 1]
}} });
</script>
<script type="text/javascript"
src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML">
</script>
<meta http-equiv="X-UA-Compatible" CONTENT="IE=EmulateIE7" />


以下テスト.
$t$を時刻,$x(t)$を状態変数,$u(t)$を入力,$y(t)$を出力として,線形のダイナミクスを
\begin{align}
\dot{x}(t) &= Ax(t) + Bu(t)\\
y(t) &= Cx(t)
\end{align}
のように記述する.
\begin{align}
A &=
\begin{bmatrix}
a_{11} & a_{12}\\
a_{21} & a_{22}
\end{bmatrix}\\
b &= \begin{bmatrix}
b_1 & b_2
\end{bmatrix}^\top
\end{align}

太字$\bm{x}$がいまいちだけど,まぁ,許容範囲かなぁ...

ふたたび再利用

よーし,また再利用しちゃうぞー

創造力

文章を「書ける人」と「書けない人」のちがい - デマこいてんじゃねえ!
高度に発達した編集は創作と区別がつかない - 雪見、月見、花見。
を読んで.

文章を書くという仕事は、ゼロを1にする作業だと思われがちだ。
(中略)
しかし実際には、文章を書くというのは100を1にする作業だ。



モノを創造するときには,インプットが大切だ.
これは多くの人が(少なくとも無意識的には)感じているはずだと思う.

そして,世の中の意見や出来事に影響されていない人は(おそらく)いないので,
厳密な意味でゼロからものを作り出すということは不可能である.

そのなかで創造とは何か?

私の意見では,100の情報を用いて何かモノを作ったときに,
100の情報が張る空間に入らない部分が十分に多いとき,それが創造である.


いいかえれば,「100を1に」したあとの1の中に,
他人では100から直感的に導けない価値が存在するとき,それが創造である.

ただし,それが世間から「創造的」だと認識されるためには高いハードルが存在する.
世間には100よりも多くの情報を持った人がいるからである.
その他人は,自分の持っている情報の張るる空間を用いて創造性を評価することになる.
作り出された価値がその空間に含まれる場合には,「何の創造性もないじゃないか」となってしまう.

このように考えると,インプットの量の重要性は明確である.
小さい空間をもとにした「創造」は世間の前では創造とみなされない.
その意味で,100を1にするという姿勢は重要な意味をもつ.

一方で,一生懸命100を1に縮めただけでは(多くの場合には)価値は生まれない.
既存の情報から新しい次元へと空間を広げる力が必要である.
その意味で,ゼロを1にする能力が必要になってくる.

したがって,この2つの能力は対立するものでなく,両方が必要なのだと私は思う.


って,考えて書いたあげく,当たり前のことしか言ってないなぁ・・・.

ふと思い立って

ふと思い立ってブログでもしてみようかと思った.
めんどくさいから,昔使ってたアカウントを再利用.
プロフィール

kawa0616

Author:kawa0616
某企業で研究開発してます.専門は制御工学と機械学習を少々.

最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
検索フォーム
RSSリンクの表示
リンク
QRコード
QR
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。