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

大まかな話は前回の解説を見てほしい。
ここでは、前回ノイズが乗っている３Dのインプットに対して、ノイズが消えるような記載をしていたが、これはそもそも入力にノイズがのっておらず、コードが間違っていたのでそこをやり直して報告する。
また、今回はノイズの大きさを変えて、フィッティング精度の違いを報告する。

やったこと

・バグフィックスしたコードの解説
・noiseの大きさを変えた時の精度について
・Gifアニメーションにしてみると、。。

・バグフィックスしたコードの解説

必要なコードは全て示すが、バグがあった３Dの部分を中心に解説する。
まず、バグがあったのは、以下のとおりである。
本来ノイズがのっている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にノイズを加えているので、以下のｙのままで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()

出力例；全体に滑らかです

以下検証です。
相関係数は関数の曲線とノイズの乗ったｙととっています。

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
ノイズが除去されているのがわかります。

同様に、平均二乗誤差の平方根も関数の曲線及びノイズの乗ったｙとの誤差を計算します。

# 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の大きさを変えた時の精度について

このままだと、ノイズの効果が今一つ不明確なので、ノイズの大きさを変えながら出力の変化を見ることとする。
最初はノイズ０の元関数

計算例.

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

・現実の系に適用しよう