大まかな話は前回の解説を見てほしい。
ここでは、前回ノイズが乗っている3Dのインプットに対して、ノイズが消えるような記載をしていたが、これはそもそも入力にノイズがのっておらず、コードが間違っていたのでそこをやり直して報告する。
また、今回はノイズの大きさを変えて、フィッティング精度の違いを報告する。
やったこと
・バグフィックスしたコードの解説
・noiseの大きさを変えた時の精度について
・Gifアニメーションにしてみると、。。
・バグフィックスしたコードの解説
必要なコードは全て示すが、バグがあった3Dの部分を中心に解説する。
まず、バグがあったのは、以下のとおりである。
本来ノイズがのっているy1でフィッティングすべきところをノイズがのっていないyでフィッティングしていた。その結果、最後の出力もノイズが全くない綺麗なカーブを描くこととなっている。
y_rbf = svr_rbf.fit(X, y).predict(X)
そこで、以下のように本来のフィッティング対象である、y1として動かしてみると、エラーが出て動かない。。。
y_rbf = svr_rbf.fit(X, y1).predict(X)
ValueError: bad input shape (100, 100)
せめて、以下のように相関などの検証をノイズがのっているオリジナルとやってみようとすると、これもエラーて動きません
rbf_corr1 = np.corrcoef(y1, test_rbf)[0, 1]
ValueError: all the input array dimensions except for the concatenation axis must match exactly
ということで、ノイズを載せる部分から変更することとして、以下のコードにたどり着きました。
最初の部分は変更ありません。
import numpy as np
from sklearn.svm import SVR
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
m = 2 #dimension
mean = np.zeros(m)
sigma = np.eye(m)
N = 100
x1 = np.linspace(-5, 5, N)
x2 = np.linspace(-5, 5, N)
X1, X2 = np.meshgrid(x1, x2)
X = np.c_[np.ravel(X1), np.ravel(X2)]
ノイズ載せる部分は、以下のように変更します。np.random.rand(X)などでごまかそうとしましたが、どうもyとの対応が悪くできません。
ということで、引数なしで乱数を与え、以下のようにy[i]に一つずつ直接加えることとしました。
ここでy.copy()はfloatなので使えないとでます。
※ここは使っていないのでこれ以上追求せず、削除しました
AttributeError: 'float' object has no attribute 'copy'
# アウトプットを算出
y = np.sin(X1).ravel()+ np.cos(X2).ravel()
for i in range(len(y)):
y[i] = y[i]+1*(0.5 - np.random.rand())
# y_o = y.copy()
y1 = y.reshape(X1.shape)
# ノイズを加える
# y1[::5] += 1 * (0.5 - np.random.rand(20,100))
以下描画部分ですが、変更なしで描画できました。
fig=plt.figure(figsize=(12, 10))
ax1 = fig.add_subplot(211, projection='3d')
ax2 = fig.add_subplot(212, projection='3d')
surf = ax2.plot_surface(X1, X2, y1, cmap='bwr', linewidth=0)
fig.colorbar(surf)
ax1.scatter3D(np.ravel(X1), np.ravel(X2), y)
ax2.set_title("Surface Plot")
ax1.set_title("scatter Plot")
plt.savefig('Surface_plot2_input.jpg')
plt.pause(1)
plt.close()
出力例
フィッティング対象も今回は直接yにノイズを加えているので、以下のyのままでOKです。
テストデータも特に変更なしです。
# フィッティング
svr_rbf = SVR(kernel='rbf', C=1e3, gamma=0.1)
y_rbf = svr_rbf.fit(X, y).predict(X)
# テストデータも準備
N = 100
test_x1 = np.linspace(-5, 5, N)
test_x2 = np.linspace(-5, 5, N)
test_X1, test_X2 = np.meshgrid(test_x1, test_x2)
test_X = np.c_[np.ravel(test_X1), np.ravel(test_X2)]
いよいよ以下で推定します。
上と同じようにsurface描画のためにtest_rbf1に変換しています。
# テストデータを突っ込んで推定してみる
test_rbf = svr_rbf.predict(test_X)
test_rbf1 = test_rbf.reshape(test_X1.shape)
fig=plt.figure(figsize=(12, 10))
ax1 = fig.add_subplot(211, projection='3d')
ax2 = fig.add_subplot(212, projection='3d')
ax1.set_title('Graph')
ax1.scatter3D(np.ravel(test_X1), np.ravel(test_X2), test_rbf)
surf = ax2.plot_surface(test_X1, test_X2, test_rbf1, cmap='bwr', linewidth=0) #Y_plot
fig.colorbar(surf)
plt.savefig('Surface_svr_rbf_prediction.jpg')
plt.pause(1)
plt.close()
以下検証です。
相関係数は関数の曲線とノイズの乗ったyととっています。
from sklearn.metrics import mean_squared_error
from math import sqrt
test_y = np.sin(test_X1).ravel() + np.cos(test_X2).ravel()
# 相関係数計算
rbf_corr = np.corrcoef(test_y, test_rbf)[0, 1]
rbf_corr1 = np.corrcoef(y, test_rbf)[0, 1]
計算例
関数との相関 Corr 0.999166
ノイズがのった入力との相関 Corr 0.960265
ノイズが除去されているのがわかります。
同様に、平均二乗誤差の平方根も関数の曲線及びノイズの乗ったyとの誤差を計算します。
# RMSEを計算
rbf_rmse = sqrt(mean_squared_error(test_y, test_rbf))
rbf_rmse1 = sqrt(mean_squared_error(y, test_rbf))
print( "RBF: RMSE %f \t\t Corr %f" % (rbf_rmse, rbf_corr))
print( "RBF: RMSE with original %f \t\t Corr %f" % (rbf_rmse1, rbf_corr1))
計算例は、それぞれ以下のとおりです。
>python svr_multi.py
RBF: RMSE 0.040633 Corr 0.999166
RBF: RMSE with original 0.286678 Corr 0.960265
両者を並べるとノイズが除去されてフィッティングされているのがわかります。
・noiseの大きさを変えた時の精度について
このままだと、ノイズの効果が今一つ不明確なので、ノイズの大きさを変えながら出力の変化を見ることとする。
最初はノイズ0の元関数
y[i] = y[i]+0*(0.5 - np.random.rand())のとき
RBF: RMSE 0.065823 Corr 0.998981
RBF: RMSE with original 0.065823 Corr 0.998981
y[i] = y[i]+0.2*(0.5 - np.random.rand())のとき
RBF: RMSE 0.002231 Corr 0.999997
RBF: RMSE with original 0.057837 Corr 0.998270
y[i] = y[i]+0.5*(0.5 - np.random.rand())のとき
RBF: RMSE 0.018703 Corr 0.999819
RBF: RMSE with original 0.142933 Corr 0.989562
y[i] = y[i]+2*(0.5 - np.random.rand())のとき
RBF: RMSE 0.091504 Corr 0.995789
RBF: RMSE with original 0.573193 Corr 0.863506
y[i] = y[i]+5*(0.5 - np.random.rand())のとき
RBF: RMSE 0.217719 Corr 0.976292
RBF: RMSE with original 1.436972 Corr 0.562649
y[i] = y[i]+10*(0.5 - np.random.rand())のとき
RBF: RMSE 0.461840 Corr 0.919116
RBF: RMSE with original 2.892038 Corr 0.341330
y[i] = y[i]+20*(0.5 - np.random.rand())のとき
RBF: RMSE 0.653215 Corr 0.837555
RBF: RMSE with original 5.740952 Corr 0.182557
y[i] = y[i]+50*(0.5 - np.random.rand())のとき
RBF: RMSE 2.042576 Corr 0.410270
RBF: RMSE with original 14.370952 Corr 0.100893
Gifアニメーションにしてみると、。。
以下のコードでGifアニメーションにしてみました。
from PIL import Image, ImageDraw
s=9
list=[0.2,0.5,1,2,5,10,20,50]
images = []
for i in list:
im = Image.open('./gifAnime/Surface_svr_rbf_prediction' + str(i)+'.jpg')
images.append(im)
images[0].save('Surface_svr_rbf_prediction.gif', save_all=True, append_images=images[1:s], duration=1000*1, loop=0)
※入力がドラスティックに変わる中、出力はかなり追従して元関数を再現しました。
入力;
まとめ
・バグフィックスしてノイズがのった入力に対してSVRの多変量回帰ができました
・ノイズを変化させたときの元関数の再現性をみました
・前回のSVRの理論からも推測できるようにノイズリダクションには有効なことが分かりました
・現実の系に適用しよう