RVMで回帰
5月も残り僅か...
新歓や講義も落ち着き始めたのでPRMLを再開。今回は7章後半の関連ベクトルマシン(Relevance Vector Machine)。
感想とか
- p57に「モデル(7.78)を用いた場合は, 計画行列Φは,...」のところで、自分はN×(N+1)の行列(1列目は全て1、残りのN×NはK_nm=k(x_n, x_{m-1})な感じの行列を想像したのだが、本文を見ると(N+1)×(N+1)の対象カーネル行列と書いてある。何か勘違いしてるのかな...？
- SVNとあんまし関係無い気g(ry
- 本題とは外れるが(C.7)のWoodBuryの公式の便利さをとようやく理解...
- p64の「各繰り返しにかかる計算時間はO(M^3)」ってのは、αiが無限大だとΣのi行i列はほぼ0になるので、モデルに含まれる基底ベクトルだけでΦを構成して、M×MのΣにしても計算結果あんま変わらなくて大丈夫→やった計算減るよ！(特にQとかSの) …みたいな感じで理解は合ってるのかな
- 超パラメータは勝手に決めてくれるけど、カーネル関数のパラメータ(ガウスカーネルとか)とかはどうすんだろ。…交差検定？
数式はちゃんと紙とペンで導出したし満足…、と言いつつほんとに動くんかいな。てなわけで以下今回書いたコード。
行列の計算が多かったので、scipyが気持ち良い。分類問題はまた明日、というか確率と画像処理論の課題ががががが。
まずは(7.87)でαの更新を行なうパターン。
αの上限値を定めたのは、値がInfになると逆行列が求められないから。
カーネル関数はガウスカーネルを用いた。
#!/usr/bin/python # -*- coding: utf-8 -*- import scipy as sp import scipy.linalg as spla import itertools as it import functools as fn class RVM(object): u""" relevance vector machine PRML 7.2.1を参考にした """ def __init__(self, kernel=lambda x, y: sp.exp(-sp.square(spla.norm(x-y))/0.05)): self._kernel = kernel def learn(self, X, t, tol=0.01, amax=1e10): u"""学習""" N = X.shape[0] a = sp.ones(N+1) # hyperparameter b = 1.0 phi = sp.ones((N, N+1)) # design matrix phi[:,1:] = [[self._kernel(xi, xj) for xj in X] for xi in X] diff = 1 while diff >= tol: sigma = spla.inv(sp.diag(a) + b * sp.dot(phi.T, phi)) m = b * sp.dot(sigma, sp.dot(phi.T, t)) gamma = sp.ones(N+1) - a * sigma.diagonal() anew = gamma / (m * m) bnew = (N - gamma.sum()) / sp.square(spla.norm(t - sp.dot(phi, m))) anew[anew >= amax] = amax adiff, bdiff = anew - a, bnew - b diff = (adiff * adiff).sum() + bdiff * bdiff a, b = anew, bnew self._a = a self._b = b self._X = X self._m = m self._sigma = sigma self._amax = amax def mean(self, x): u"""予測値の平均""" ret = 0 phi = sp.append([1.0], [self._kernel(x, xi) for xi in X]) for i in range(len(self._m)): if self._a[i] < self._amax: ret += self._m[i] * phi[i] return ret def variance(self, x): u"""予測値の分散""" phi = sp.append([1.0], [self._kernel(x, xi) for xi in X]) return 1.0/self._b + sp.dot(phi.T, sp.dot(self._sigma, phi)) def _get_a(self): return self._a # 超パラメータ a = property(_get_a) def _get_rv_index(self): return self._a[1:] < self._amax # 関連ベクトルのインデックス rv_index = property(_get_rv_index) if __name__ == '__main__': # めんどうなのでデータはここに;-p data = sp.array([ [0.000000, 0.349486], [0.111111, 0.830839], [0.222222, 1.007332], [0.333333, 0.971507], [0.444444, 0.133066], [0.555556, 0.166823], [0.666667, -0.848307], [0.777778, -0.445686], [0.888889, -0.563567], [1.000000, 0.261502], ]) X = data[:, 0] t = data[:, 1] rvm = RVM() rvm.learn(X, t) print "a=", rvm.a # 描画 import matplotlib.pyplot as plt x = sp.linspace(0, 1, 50) # +-標準偏差1つ分の幅の表示 meshx, meshy = sp.meshgrid(sp.linspace(0, 1, 200), sp.linspace(-1.5, 1.5, 200)) meshz = [[abs(rvm.mean(meshx[j][i]) - meshy[j][i]) <= sp.sqrt(rvm.variance(meshx[j][i])) for i in range(len(meshx[0]))] for j in range(len(meshx))] plt.contour(meshx, meshy, meshz, 1) plt.spring() # 入力 plt.scatter(X, t, label="input") # 関連ベクトルの描画 plt.scatter(X[rvm.rv_index], t[rvm.rv_index], marker='d', color='r', label="relevance vector") # 元の曲線を表示 y = sp.sin(2*sp.pi*x) plt.plot(x, y, label="sin(2pix)") # RVMの予測の平均 y_ = [rvm.mean(xi) for xi in x] plt.plot(x, y_, label="RVM regression") # 表示範囲の調整 plt.xlim(-0.05, 1.05) plt.ylim(-1.5, 1.5) # label表示 plt.legend() plt.show()
出力は
a= [ 1.00000000e+10 1.00000000e+10 1.00000000e+10 8.99348648e-01 1.00000000e+10 1.00000000e+10 1.00000000e+10 1.00000000e+10 1.95794077e+00 1.00000000e+10 1.09871526e+01]
と確かに疎。
グラフは以下。ピンクの線は平均から標準偏差1つ分に相当する幅で等高線で引いたもの。赤い点は関連ベクトル。
…上手く動いてるようなので、逐次的疎ベイジアン学習アルゴリズムとやらを試してみる。
#!/usr/bin/python # -*- coding: utf-8 -*- import scipy as sp import scipy.linalg as spla import itertools as it import functools as fn class RVM(object): u""" relevance vector machine PRML 7.2.2を参考にした """ def __init__(self, kernel=lambda x, y: sp.exp(-sp.square(spla.norm(x-y))/0.05)): self._kernel = kernel def learn(self, X, t, tol=0.01, amax=1e10): u"""学習""" N = X.shape[0] a = sp.ones(N+1)*amax # hyperparameter phiT = sp.ones((N+1, N)) # design matrix(*transposed*) phiT[1:] = [[self._kernel(xi, xj) for xj in X] for xi in X] bases = set([0]) # basis func indexes a[0] = 1.0 b = 1 diff = 1 while diff >= tol: indexes = list(bases) phiT_ = phiT[indexes] sigma = spla.inv(sp.diag(a[indexes]) + b * sp.dot(phiT_, phiT_.T)) m = b * sp.dot(sigma, sp.dot(phiT_, t)) S = Q = b * phiT - b*b*sp.dot(phiT, sp.dot(phiT_.T, sp.dot(sigma, phiT_))) Q = sp.dot(Q, t) S = sp.dot(S, phiT.T).diagonal() q = Q / (1 - S / a) s = S / (1 - S / a) anew = a.copy() for i in range(N+1): q2 = q[i]*q[i] if q2 > s[i] and a[i] < amax: anew[i] = s[i]*s[i]/(q2-s[i]) elif q2 > s[i] and a[i] >= amax: bases.add(i) anew[i] = s[i]*s[i]/(q2-s[i]) elif q2 <= s[i] and a[i] < amax: bases.remove(i) anew[i] = amax bnew = (sp.dot(a[indexes], sigma.diagonal())+N-len(indexes)) / sp.square(spla.norm(t - sp.dot(phiT_.T, m))) adiff, bdiff = anew - a, bnew - b diff = (adiff * adiff).sum() + bdiff * bdiff a, b = anew, bnew self._a = a self._b = b self._m = m self._sigma = sigma self._indexes = list(bases) def mean(self, x): u"""予測値の平均""" phi = sp.append([1.0], [self._kernel(x, xi) for xi in X])[self._indexes] return sp.dot(phi, self._m) def variance(self, x): u"""予測値の分散""" phi = sp.append([1.0], [self._kernel(x, xi) for xi in X])[self._indexes] return 1.0/self._b + sp.dot(phi.T, sp.dot(self._sigma, phi)) def _get_a(self): return self._a a = property(_get_a) def _get_b(self): return self._b_ b = property(_get_b) def _get_m(self): return self._a m = property(_get_m) def _get_rv_index(self): indexes = [x-1 for x in self._indexes] if len(indexes) > 0 and indexes[0] == -1: indexes = indexes[1:] return indexes # 関連ベクトルのインデックス rv_index = property(_get_rv_index) if __name__ == '__main__': # めんどうなのでデータはここに;-p data = sp.array([ [0.000000, 0.349486], [0.111111, 0.830839], [0.222222, 1.007332], [0.333333, 0.971507], [0.444444, 0.133066], [0.555556, 0.166823], [0.666667, -0.848307], [0.777778, -0.445686], [0.888889, -0.563567], [1.000000, 0.261502], ]) X = data[:, 0] t = data[:, 1] rvm = RVM() rvm.learn(X, t) print "a=", rvm.a # 描画 import matplotlib.pyplot as plt x = sp.linspace(0, 1, 50) # +-標準偏差1つ分の幅の表示 meshx, meshy = sp.meshgrid(sp.linspace(0, 1, 200), sp.linspace(-1.5, 1.5, 200)) meshz = [[abs(rvm.mean(meshx[j][i]) - meshy[j][i]) <= sp.sqrt(rvm.variance(meshx[j][i])) for i in range(len(meshx[0]))] for j in range(len(meshx))] plt.contour(meshx, meshy, meshz, 1) plt.spring() # 入力 plt.scatter(X, t, label="input") # 関連ベクトルの描画 plt.scatter(X[rvm.rv_index], t[rvm.rv_index], marker='d', color='r', label="relevance vector") # 元の曲線を表示 y = sp.sin(2*sp.pi*x) plt.plot(x, y, label="sin(2pix)") # RVMの予測の平均 y_ = [rvm.mean(xi) for xi in x] plt.plot(x, y_, label="RVM regression") # 表示範囲の調整 plt.xlim(-0.05, 1.05) plt.ylim(-1.5, 1.5) # label表示 plt.legend() plt.show()
出力は以下。先程と殆どかわらず。ほぇー。
a= [ 1.00000000e+10 1.00000000e+10 1.00000000e+10 8.99430644e-01 1.00000000e+10 1.00000000e+10 1.00000000e+10 1.00000000e+10 1.96108669e+00 1.00000000e+10 1.10242579e+01]
グラフは
ちなみにP59に「データが存在しない領域で予測分散が小さく...」と書いてあったので、データを一部削除して試してみた。
…確かに図6.8のガウス過程さんに比べて残念っぽい。
