Edited at

【機械学習入門】「回帰分析;SVRとかとか…二次元から多変量回帰」で遊んでみた♬ノイズリダクションに有効

大まかな話は前回の解説を見てほしい。

ここでは、前回ノイズが乗っている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の理論からも推測できるようにノイズリダクションには有効なことが分かりました

・現実の系に適用しよう