Zenn
🍣

Gpyoptによるベイズ最適化でフィッティング

2022/07/22に公開
Python
gpyopt
tech

GPyOptによるベイズ最適化＋実験計画 (1) ~2変数~

参考にしたリンク

ライブラリ

  • python 3.8.13
  • matplotlib 3.1.3
  • gpy 1.10.0
  • gpyopt 1.2.6
  • numpy 1.22.3
  • pandas 1.2.5

(注) GPyOptのplot_aquisition()でエラーが発生。Matplotlib 3.1.3以前では回避できるとのことから各ライブラリをダウングレードして使用。

やりたいこと

  • ax^2+bの放物線プロットを作成(模擬実験データ)
  • ax^2+bと仮定してある(a,b)の組に対し計算上のyとデータのyの差の二乗を計算。これを目的変数とする。
  • ベイズ最適から(a,b)の組に対し、最も目的変数の小さい組み合わせを算出する。すなわち(a, b)のフィティングを行う

コード

# 必要なライブラリをインポート
import GPy
import GPyOpt
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline

# 元となる関数を定義
def func_test(a,b,x):
  y = a*x**2+b
  return y

# a=2.0, b = 1.3として、xとyの関係を求める。
a= 2.0
b= 1.3
# 下記x_expとy_expを実験データとして扱う。
x_exp = np.linspace(-1.0,1.0,100)
y_exp = func_test(a,b,x_exp)
plt.scatter(x_exp, y_exp)



# 誤差を算出する関数を作る
def func_def(a,b,x_exp,y_exp):
  y_calc = func_test(a,b,x_exp) #x_expがndarrayなのでここもndarrayになる。
  delta_sq = (y_calc - y_exp)**2
  sum_delta = np.sum(delta_sq)
  return sum_delta

# 初点を(a,b)=(0,0)とする。2D arrayであることに注意。
ab_pred = np.asarray([[0.0, 0.0]])
# (a,b)=(0,0)のときのδ^2を算出する。2D arrayのデータを取り出して計算後、再度2D arrayとして保存
delta_calc = np.asarray([[func_def(ab_pred[0,0],ab_pred[0,1],x_exp,y_exp)]])

まずは、初期データから次の予測ができることを確認。

# 保有データから次のデータの予測を1回　行う
bounds = [{'name': 'a', 'type': 'continuous', 'domain': (-5.0, 5.0)},
          {'name': 'b', 'type': 'continuous', 'domain': (-5.0, 5.0)}]

params = {'acquisition_type':'EI',#獲得関数としてEIを指定    
          'f':None, #最適化する関数の設定（実験結果は分からないので設定しない）．
          'domain':bounds, #パラメータの探索範囲の指定
          'model_type':'GP', #標準的なガウシアンプロセスを指定．
          'X':ab_pred, #既知データの説明変数（インプットするx）
          'Y':delta_calc, #既知データの目的変数（インプットするy）
          'de_duplication':True, #重複したデータをサンプルしないように設定．
          "normalize_Y":True, #defaltはTrue
          "exact_feval":False #defaltはFalse
          } 
myBopt = GPyOpt.methods.BayesianOptimization(**params)
myBopt.plot_acquisition()

x_suggest = myBopt.suggest_next_locations(ignored_X=ab_pred) 
y_predict = myBopt.model.model.predict(x_suggest) #y_predictは(予測平均，予測分散)がタプルで返ってくる
y_mean=y_predict[0]
y_variance=y_predict[1]

print(f'つぎは{x_suggest}、yの予測平均は{y_mean}、分散は{y_variance}')


つぎは[[-3.02910399  1.09674807]]、yの予測平均は[[0.]]、分散は[[1.e-09]]

と出ればOK。

# 更新用の2d array作成。初期データに推奨a,bのペア、および新規目的変数を作成
ab_new_array = np.vstack((ab_pred,x_suggest)) # np.stack(())が二重になる点注意。
delta_new = np.asarray([[func_def(x_suggest[0,0],x_suggest[0,1],x_exp,y_exp)]])
delta_new_array = np.vstack((delta_calc,delta_new)) # np.stack(())が二重になる点注意。

# 更新データから次のデータの予測を行う
# 新しいパラメータ
num = 10 #追加回数
for n in range(num):
  params_2 = {'acquisition_type':'EI',#獲得関数としてEIを指定    
            'f':None, #最適化する関数の設定（実験結果は分からないので設定しない）．
            'domain':bounds, #パラメータの探索範囲の指定
            'model_type':'GP', #標準的なガウシアンプロセスを指定．
            'X':ab_new_array, #既知データの説明変数（インプットするx）ここが変更
            'Y':delta_new_array, #既知データの目的変数（インプットするy）ここが変更
            'de_duplication':True, #重複したデータをサンプルしないように設定．
            "normalize_Y":True, #defaltはTrue
            "exact_feval":False, #defaltはFalse
            'maximize':False
            } 
  myBopt_2 = GPyOpt.methods.BayesianOptimization(**params_2)
  ab_suggest_2 = myBopt_2.suggest_next_locations(ignored_X=ab_new_array) 
  delta_predict_2 = myBopt_2.model.model.predict(ab_suggest_2) 
  delta_mean_2 =delta_predict_2[0]
  delta_variance_2=delta_predict_2[1]
  ab_new_array = np.vstack((ab_new_array, ab_suggest_2)) # np.stack(())が二重になる点注意。
  delta_new = np.asarray([[func_def(ab_suggest_2[0,0],ab_suggest_2[0,1],x_exp,y_exp)]])
  delta_new_array = np.vstack((delta_new_array,delta_new)) # np.stack(())が二重になる点注意。
  
  print(f'つぎは{ab_suggest_2}、deltaの予測平均は{delta_mean_2}、分散は{delta_variance_2},実際は{delta_new}')


myBopt_2.plot_acquisition()

Discussion