本記事では、目的関数Yの最小値または最大値と、説明変数Xを求めるための雛形コードを載せました。ベイズ最適化のPythonライブラリは複数あります。本記事ではGPyOptについてで、インストールはcondaで次のようにします。
conda install -c conda-forge gpyopt
はじめに、最適化する出力指標を目的関数と呼び、入力因子のことを説明変数と言います。例えば、下図のようなz=x(exp(-x/5)^2-(y/5)^2)の場合は、zが目的関数でxとyが説明変数です。このzの最大値を求めたいとします。期待する結果は(x, y)=(3.5前後, 0)で最大値は2前後です。(参考までに説明変数は3つ以上でも出来ます。また制約条件を設けることも可能です。本コード中にコメントアウト文で例を記載しています)
本コードを実行すると、(x, y)の探索履歴を下図のようにグラフで出力します。探索回数が増えるにつれて、〇の色の濃淡が濃くなるように表示しています。下図3つの内の一番上図は、9回探索した結果で(x, y)=(3.5, 0)近傍を多く探索している様子です。その次の中央の図は、19回探索した結果で(x, y)=(3.5, 0)付近にプロットがあり最大値を見つけたようです。一番下の図は、一応念のため他の箇所(x, yの各々の探索範囲の上下限)も探索している様子です。ここまですると普通にDOEしているのと変わらないような状況です。
そして、下図左のように最適結果、下図右のように探索履歴をcsvファイルに保存します。
参考までに、Jupyter Lab等での使用環境の場合、探索結果は下図のようにたった一行で可視化できます。
■本プログラム
本コードは、Jupyter labで作成後に.pyファイルに出力して一部編集したものです。Anaconda PromptのようなCUIで実行できます。
import numpy as np
import GPy
import GPyOpt
def f(x):
x0, x1 = x[:,0], x[:,1]
return x0*(np.exp(-1 * pow(x0/5, 2) - pow(x1/5, 2)))
bounds = [
{'name': 'x0', 'type': 'continuous', 'domain': (-10, 10)},
{'name': 'x1', 'type': 'discrete', 'domain': (-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10)}
]
'''
# 制約条件を設定したい場合
constraints = [
{
"name": "constr_1",
"constraint": "(x[:,0] + x[:,1]) - 10 - 0.001" # x0 + x1 -10 <= 0.001 の制約設定したい場合
},
{
"name": "constr_2",
"constraint": "10 - (x[:,0] + x[:,1]) - 0.001" # x0 + x1 -10 >= -0.001 の制約設定したい場合
}
]
'''
flag_maximize = True
num = 30
myBopt = GPyOpt.methods.BayesianOptimization(
f=f,
domain=bounds,
maximize = flag_maximize,
model_type = 'GP'
)
myBopt.run_optimization(max_iter = num)
if flag_maximize:
opt_fx = myBopt.fx_opt * (-1)
else:
opt_fx = myBopt.fx_opt
opt_x = myBopt.x_opt
print(opt_fx)
print(opt_x)
columns_name_list = ['x0', 'x1']
columns_name_list_str = ','.join(columns_name_list)
columns_name_list_str = columns_name_list_str + ',Y' + '\n'
with open('opt_data.csv', 'w') as f:
f.write(columns_name_list_str)
buf_list = []
for X in opt_x.tolist():
buf_list.append(X)
buf_list.append(opt_fx)
opt_list = list(map(str, buf_list))
opt_str = ','.join(opt_list)
f.write(opt_str + '\n')
import matplotlib.pyplot as plt
import japanize_matplotlib
step = int(num/3)
for i in range(step-1, num, step):
history_X = myBopt.X[0:i]
plt.rcParams['font.size'] = 13
plt.figure()
plt.title('開始から' + str(i) + '番目までの探索履歴')
plt.xlabel("X")
plt.ylabel("Y")
plt.xlim(-11, 11)
plt.ylim(-11, 11)
sc = plt.scatter(
history_X[:,0],
history_X[:,1],
s = 100,
c = range(len(history_X)),
cmap = "Greys")
plt.colorbar(sc)
plt.plot(history_X[:,0], history_X[:,1], linestyle=":", color='k')
plt.grid()
plt.savefig('history_' + str(i) + '.png')
plt.close()
if flag_maximize:
history_Y = myBopt.Y[0:num] * (-1)
else:
history_Y = myBopt.Y[0:num]
history_Y_list = history_Y.tolist()
history_X_list = history_X.tolist()
with open('history_data.csv', 'w') as f:
f.write(columns_name_list_str)
for X, Y in zip(history_X_list, history_Y_list):
X.extend(Y)
buf_list = list(map(str, X))
buf_str = ','.join(buf_list)
f.write(buf_str + '\n')
以上
<広告>
リンク