가장 쉬운 베이지안 옵티마이제이션 [Python] by 바죠

가장 쉬운 베이지안 옵티마이제이션(Bayesian optimization)

H: 가설
D: 데이터

베이즈 추정의 정의는 사전확률을 행동의 관찰(정보)에 의거해 사후확률로 베이즈 업데이트 하는 것이다.
추론 대상의 사전 확률과 추가적인 정보를 통해 해당 대상의 사후 확률을 추론하는 방법이다. 
베이즈 추정에는 축차합리성이 있다.

첫번째 정보와 두번째 정보를 동시에 고려하는 것 == 첫번째 정보를 고려하여 사후확률을 계산한 후, 이를 두번째 정보에 대한 사전확률로 사용하는 것 == 학습이 진행된다고 볼 수 있다.


베이즈 확률론(Bayesian probability)은 확률을 '지식 또는 믿음의 정도를 나타내는 양'으로 해석하는 확률론이다.
확률을 발생 빈도(frequency)나 어떤 시스템의 고유한 물리적 속성으로 여기는 것과는 다른 해석이다. 
18세기 통계학자 토머스 베이즈의 이름에서 왔다.
어떤 가설의 확률을 평가하기 위해서 사전 확률을 먼저 밝히고 새로운 관련 데이터에 의한 새로운 확률값을 변경한다.
베이즈 통계학 (Bayesian statistics)은 하나의 사건에서의 믿음의 정도 (degree of belief)를 확률로 나타내는 베이즈 확률론에 기반한 통계학 이론이다. 
믿음의 정도는 이전 실험에 대한 결과, 또는 그 사건에 대한 개인적 믿음 등, 그 사건에 대한 사전 지식에 기반할 수 있다. 
이것은 확률을 많은 시도 후의 사건의 상대적 빈도의 극한으로 보는 빈도주의자 (frequentist) 등 많은 다른 확률에 대한 해석과는 차이가 있다. 
어떠한 경우는 동일한 해석을 주기도 한다.

cos(x)라는 함수를 생각해 보자. 
여기서 우리는 cos(x)를 예로 들었다. 
우리의 가정하나를 잊지 말아야 한다. 이 함수 계산이 엄청나게 시간이 많이 걸린다고 가정하고 있다.
유한한 개수의 실측을 통해서 얻을 수 있는 함수값들을 생각한다. 
당연히 함수꼴에 대한 한정적인 정보만을 취할 수밖에 없다.
적절히 분포되어 있는 실측 함수값들이 있다고 가정한다.
4, 8, 16 points.



기계학습에서 사용되는 모델의 성능 측정에서 cross validation score 함수는 연속함수일리가 없다. 이 함수 계산은 제법 시간이 걸리는 계산이다. 데이터가 많을 수록 더욱더 오랜 시간 계산을 수행할 것이다.
유한한 데이터 기반의 정밀도 측정은 불연속 함수를 사실상 가정하고 있다.
함수의 미분은 해석적으로 얻기 어렵다. 당연하다. 

노이즈를 포함한 함수값이 얻어진 경우들, 즉, 데이터를 확보한 경우들에 대한 철저한 고찰이 필요하다. 
계속해서 모아둔다. 이들 데이터를 활용한다. 
데이터로부터 더 최적화된 함수값을 예측한다.
점진적으로 데이터를 확충해나간다.  축차합리성을 이용한다.



베이지안 옵티마이제이션(Bayesian optimization)
매우 복잡한 문제에 적용할 수 있는 아주 일반적인 최적화 알고리듬은 아니다.
하지만, 적당한 크기의 문제에 대해서 적절한 최적화 방법이라고 볼 수 있다.
특히, 기계학습의 경우, hyperparameter 최적화에 활용되고 있다.
이 문제의 크기(hyperparameter들의 수)와 베이지안 옵티마이제이션는 궁합이 잘 맞는 경우라고 볼 수 있다.
이 경우, 기계학습 문제(hyperparameter 결정 문제)를 기계학습 방법(베이지안 옵티마이제이션, Bayesian optimization)으로 풀어내는 것이다. 
최적화해야 할 변수의 숫자가 20개 이하라고 볼 수 있다. 

베이지안 옵티마이제이션(Bayesian optimization) 방법의 일반적 성질들:
1. 순차적 접근법이다. 계산이 병렬화되지 않는다.

2. 목적함수의 도함수를 이용하지 않는다. [도함수를 알 수 없을 정도로 이산적인 경우에 적합하다. 
문제에 따라서는 도함수가 정의되지 않을 수 있다. 
물론 도함수가 정의되어도 계산하는 것이 복잡할 때에도 적용할 수 있다. 이 경우, 자동구배를 활용할 수도 있겠다.]

3. 기계학습 방법을 이용해서 보다 더 좋은 해가 있는 곳을 예측한다. (surrogate model; 대리 모델)[인공지능을 이용한다.]

4. 여러 가지 model이 가능하다. 또한 model 선택에 민감하게 결과가 변화할 수 있다. 

5. 목적함수 계산을 너무 많이 할 수 없는 상황에 적절한 방법으로 알려져 있다. 
[즉, 계산량이 너무 많을 경우에 해당한다. 목적함수가 아주 비싼 대가를 치르는 경우는 무수히 많다.]
다시 말해서, 기계학습 방법으로 surrogate model을 만드는 비용이 아주 싸게 먹히는 경우에 적합한 최적화 방법이다. 
(하지만, 데이터 수가 너무 많을 경우 문제가 될 수 있다. 
데이터 갯수의 삼승에 비례하는 문제가 발생한다. 이것은 역행렬 계산과 연관이 있는 복잡도이다.)
목적함수 계산에 얼마나 많은 시간이 소모되는지? 기대 향상을 얻어내는데 걸리는 시간을 직접 비교해 보아야 한다.

6. 목적함수가 노이즈를 가지고 있을 때 사용할 수 있는 일반적인 최적화 알고리듬이다.

7. 추가적인 데이터의 보강이 이루어 질 때, 지속적으로 기계학습 방법을 이용할 수 있다. 

일반적인 국소 최소화 알고리듬과는 용도가 다르다는 점에 유의해야한다.
예를 들어, 아주 간단한 함수 국소 최소화에는 여전히 전통적인 국소 최소화 알고리듬이 더 뛰어나다. [Nelder-Mead]
해석적으로 목적함수의 도함수가 알려진 경우는 전통적인 계산 방법을 사용해야한다. [BFGS]
기계학습 방법은 좀 더 상황이 꼬인 경우에 사용하는 것이다. 예를 들어, 노이즈가 있는 경우를 들 수 있다.

기대 향상을 아래와 같이 정의한다.
현재까지 최고의 함수값으로 부터 더 올라갈 수 있는 정도로 정의한다. 
최고값, 최고의 해를 알고 있다. 아울러, 기존의 데이터를 통해서 학습한 모델을 가지고 있다. 
이 모델이 예측하는 것을 기준으로 일을 진행한다. 
새로운 해를 가정하고 예측값과 분산을 모델로부터 얻어낸다.
예측값과 분산을 바탕으로 기대 향상을 얻어낸다. 
정규분포로 부터 정의되는 cdf, pdf 함수를 사용한다.
xi=0.01로 둔다.

통상 세가지 중 한 가지를 사용하여 예측을 수행한다.
maximum probability of improvement (MPI)
expected improvement (EI)
upper confidence bound (UCB)

\[ f(x+) \]:  현재까지 최고값.
\[ x+ \]  :  현재까지 최고의 해
\[\xi = 0.01\] 로 고정한다. 대리 모델을 활용한 기대 향상의 예측을 생각해 볼 수 있다. 
\[ f(\vec{x^{+}})  \]
\[ \vec{x^{+}}  \]
\[ \xi=0.01  \]

먼저 임의로 몇 개의 점들에 대해서 목적함수들을 계산한다. 

\[ \{{\vec{x}\}  } \]
\[ \{ {f(\vec{x})\} } \]

최고의 기대 향상을 얻기 위해서 또 다른 보조 최적화 방법을 활용한다.
수치적 미분을 활용하여 보조 최적화 작업을 마무리 할 수 있다. 이에 따른 계산량이 방대하지 않기 때문이다. 

\[ EI(\vec{x})   =  \mathbb{E}_{max}\{ f(\vec{x}) -f(\vec{x}^{+}),0 \}\]

---------------------------------------------------------------------------------------------------------------------
import scipy.optimize as optimize
from bayes_opt import BayesianOptimization
fun = lambda x: (x[0] - 1)**2 + (x[1] - 2.5)**2
res = optimize.minimize(fun, (2, 0), method='TNC', tol=1e-10)
print(res.x)
print(res.fun)
def fun(x1,x2):
      return -((x1-1)**2+ (x2-2.5)**2)
BO = BayesianOptimization(fun, {'x1': (0., 2.), 'x2': (0., 3.) }, verbose=2)
BO.maximize(init_points=2, n_iter=20)
for i, res in enumerate(BO.res):
    print("Iteration {}: nt{}".format(i, res))
print(BO.max)

[1.  2.5]
3.616840142797934e-17
|   iter    |  target   |    x1     |    x2     |
-------------------------------------------------
|  1        | -0.4267   |  1.445    |  2.979    |
|  2        | -2.982    |  1.061    |  0.7742   |
|  3        | -1.093    |  0.7062   |  1.497    |
|  4        | -2.886    |  1.206    |  0.8138   |
|  5        | -0.01511  |  0.8914   |  2.443    |
|  6        | -1.153    |  0.0      |  2.891    |
|  7        | -1.077    |  2.0      |  2.222    |
|  8        | -1.317    |  0.0      |  1.937    |
|  9        | -0.286    |  0.8103   |  3.0      |
|  10       | -0.1592   |  1.221    |  2.168    |
|  11       | -0.03283  |  1.158    |  2.588    |
|  12       | -0.105    |  0.9055   |  2.19     |
|  13       | -0.01146  |  1.056    |  2.409    |
|  14       | -0.0187   |  0.9762   |  2.635    |
|  15       | -0.000210 |  1.006    |  2.513    |
|  16       | -0.001221 |  0.9822   |  2.47     |
|  17       | -0.002044 |  0.9604   |  2.522    |
|  18       | -0.000427 |  1.019    |  2.509    |
|  19       | -9.666e-0 |  1.007    |  2.493    |
|  20       | -8.04e-05 |  1.008    |  2.504    |
|  21       | -0.000309 |  0.9914   |  2.485    |
|  22       | -0.000311 |  1.013    |  2.512    |
=================================================

{'target': -8.039578136899103e-05, 'params': {'x1': 1.0078913215228398, 'x2': 2.5042570912595523}}

---------------------------------------------------------------------------------------------------------------------
from bayes_opt import BayesianOptimization
import numpy as np
def target(x):
    return np.exp(-(x-3)**2) + np.exp(-(3*x-2)**2) + 1/(x**2+1)
bayes_optimizer = BayesianOptimization(target, {'x': (-2, 6)}, random_state=0)
bayes_optimizer.maximize(init_points=2, n_iter=14, acq='ei', xi=0.01)

|   iter    |  target   |     x     |
-------------------------------------
|  1        |  0.8386   |  2.391    |
|  2        |  0.6615   |  3.722    |
|  3        |  0.2      | -2.0      |
|  4        |  0.5116   |  1.215    |
|  5        |  0.02715  |  6.0      |
|  6        |  1.057    |  2.751    |
|  7        |  0.9193   | -0.2969   |
|  8        |  1.1      |  2.996    |
|  9        |  1.132    |  0.2321   |
|  10       |  1.05     |  0.1093   |
|  11       |  1.498    |  0.4547   |
|  12       |  1.707    |  0.6296   |
|  13       |  0.08147  |  4.796    |
|  14       |  0.4498   | -1.106    |
|  15       |  1.662    |  0.7032   |
|  16       |  0.03633  |  5.399    |
=====================================


---------------------------------------------------------------------------------------------------------------------
XGBoost의 hyperparameter 최적화 그리고 베이지안 옵티마이제이션
import numpy as np
from numpy import loadtxt
from xgboost import XGBClassifier
from bayes_opt import BayesianOptimization
from sklearn.model_selection import cross_val_score
pbounds = {  'learning_rate': (0.01, 1.0),  'n_estimators': (100, 1000), 'max_depth': (3,10),   
    'subsample': (1.0, 1.0),  # Change for big datasets
    'colsample': (1.0, 1.0),   # Change for datasets with lots of features
    'gamma': (0, 5)}
def xgboost_hyper_param(learning_rate, n_estimators, max_depth, subsample, colsample, gamma):
    dataset = loadtxt('pima-indians-diabetes.csv', delimiter=",")
    X = dataset[:,0:8]
    y = dataset[:,8]
    max_depth = int(max_depth)
    n_estimators = int(n_estimators)
    clf = XGBClassifier( max_depth=max_depth, learning_rate=learning_rate, n_estimators=n_estimators, subsample=subsample, colsample=colsample, gamma=gamma)
    return np.mean(cross_val_score(clf, X, y, cv=3, scoring='roc_auc'))
optimizer = BayesianOptimization( f=xgboost_hyper_param, pbounds=pbounds, random_state=1)
optimizer.maximize(init_points=3, n_iter=24, acq='ei', xi=0.01)

---------------------------------------------------------------------------------------------------------------------
import numpy as np
from numpy import loadtxt
from xgboost import XGBClassifier
from bayes_opt import BayesianOptimization
from sklearn.model_selection import cross_val_score
pbounds = {
    'learning_rate': (0.01, 1.0),
    'n_estimators': (100, 1000),
    'max_depth': (3,10),
    'subsample': (1.0, 1.0),  # Change for big datasets
    'colsample': (1.0, 1.0),  # Change for datasets with lots of features
    'gamma': (0, 5)}
def xgboost_hyper_param(learning_rate, n_estimators, max_depth, subsample, colsample, gamma):
    max_depth = int(max_depth)
    n_estimators = int(n_estimators)
    clf = XGBClassifier( max_depth=max_depth, learning_rate=learning_rate, n_estimators=n_estimators,
         subsample=subsample, colsample=colsample, gamma=gamma)
    return np.mean(cross_val_score(clf, X, y, cv=3, scoring='roc_auc'))
def boexe(X,y):
    optimizer = BayesianOptimization( f=xgboost_hyper_param, pbounds=pbounds, random_state=1)
    optimizer.maximize(init_points=3, n_iter=50, acq='ei', xi=0.01)
dataset = loadtxt('pima-indians-diabetes.csv', delimiter=",")
X = dataset[:,0:8]
y = dataset[:,8]
boexe(X,y)

---------------------------------------------------------------------------------------------------------------------
import lightgbm as lgb
from bayes_opt import BayesianOptimization
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_boston
from sklearn.metrics import accuracy_score,confusion_matrix
import numpy as np
def lgb_evaluate(numLeaves, maxDepth, scaleWeight, minChildWeight, subsample, colSam):
    reg=lgb.LGBMRegressor(num_leaves=31,
                          max_depth= 2,
                          scale_pos_weight= scaleWeight,
                          min_child_weight= minChildWeight,
                          subsample= 0.4,
                          colsample_bytree= 0.4,
                         learning_rate=0.05,
                         n_estimators=20)
#   scores = cross_val_score(reg, train_x, train_y, cv=5, scoring='roc_auc')
    scores = cross_val_score(reg, train_x, train_y, cv=5, scoring='neg_mean_squared_error')
    return np.mean(scores)
def bayesOpt(train_x, train_y):
    lgbBO = BayesianOptimization(lgb_evaluate, {                                               
                                                'numLeaves':  (5, 90),
                                                'maxDepth': (2, 90),
                                                'scaleWeight': (1, 10000),
                                                'minChildWeight': (0.01, 70),
                                                'subsample': (0.4, 1),                                               
                                                'colSam': (0.4, 1)
                                            })
    lgbBO.maximize(init_points=5, n_iter=50)
    print(lgbBO.res)
boston = load_boston()
X, y = boston.data, boston.target
train_x, X_test, train_y, y_test = train_test_split(X, y, test_size=0.2)
bayesOpt(train_x, train_y)

---------------------------------------------------------------------------------------------------------------------
import lightgbm as lgb
from bayes_opt import BayesianOptimization
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
from numpy import loadtxt
from sklearn.metrics import accuracy_score,confusion_matrix
import numpy as np
def lgb_evaluate(numLeaves, maxDepth, scaleWeight, minChildWeight, subsample, colSam):
    clf = lgb.LGBMClassifier(
        objective = 'binary',
        metric= 'auc',
        reg_alpha= 0,
        reg_lambda= 2,
#       bagging_fraction= 0.999,
        min_split_gain= 0,
        min_child_samples= 10,
        subsample_freq= 3,
#       subsample_for_bin= 50000,
#       n_estimators= 9999999,
        n_estimators= 99,
        num_leaves= int(numLeaves),
        max_depth= int(maxDepth),
        scale_pos_weight= scaleWeight,
        min_child_weight= minChildWeight,
        subsample= subsample,
        colsample_bytree= colSam,
        verbose =-1)
    scores = cross_val_score(clf, train_x, train_y, cv=5, scoring='roc_auc')
    return np.mean(scores)
def bayesOpt(train_x, train_y):
    lgbBO = BayesianOptimization(lgb_evaluate, {                                               
                                                'numLeaves':  (5, 90),
                                                'maxDepth': (2, 90),
                                                'scaleWeight': (1, 10000),
                                                'minChildWeight': (0.01, 70),
                                                'subsample': (0.4, 1),                                               
                                                'colSam': (0.4, 1)
                                            })
    lgbBO.maximize(init_points=5, n_iter=50)
    print(lgbBO.res)
dataset = loadtxt('pima-indians-diabetes.csv', delimiter=",")
X = dataset[:,0:8]
y = dataset[:,8]
train_x, X_test, train_y, y_test = train_test_split(X, y, test_size=0.2)
bayesOpt(train_x, train_y)

---------------------------------------------------------------------------------------------------------------------
scikit-learn에서 제공하는 다양한 분류, 회귀 방법들이 있다. 확장성을 보유한 의사결정 나무(decision tree)를 사용하는 방법이 독특하다. 나아가, 나무를 여러 개를 사용하는 방법이 개발되어 있다. bootstrap aggregating(bagging)이라는 방법의 하나로 random forest가 있다. 이것은 부트스트래핑 (bootstraping)을 사용한다. 즉, 중복을 허용하는 resampling을 시도한다. 
따라서, 우리는 주어진 데이터보다 더 많은 재추출된 데이터를 가지게 된다. 

다수의 의사결정 나무들을 가질수 있게 된다. 이를 숲이라고 한다. 임의로 resampling을 했기 때문에 random forest라고 이름을 붙였다. 즉, 많은 수의 sub-sample들을 활용하는 것이 특징이다. 아주 성능이 뛰어난 분류, 회귀 방법이다. 물론, LightGBM, XGBoost 등 보다 다소 성능이 좋지 못하지만 많은 연구에 활용되었다. 현재로서는 LightGBM, XGBoost, CATboost 등이 많이 활용되고 있다. 특히, LightGBM은 아주 빠른 실행 속도를 가지고 있다.

---------------------------------------------------------------------------------------------------------------------

---------------------------------------------------------------------------------------------------------------------

#!/usr/bin/env python
# coding: utf-8

get_ipython().run_line_magic('matplotlib', 'inline')


# # Tuning a scikit-learn estimator with `skopt`
# Gilles Louppe, July 2016
# Katie Malone, August 2016
# Reformatted by Holger Nahrstaedt 2020
# .. currentmodule:: skopt
# If you are looking for a :obj:`sklearn.model_selection.GridSearchCV` replacement checkout
# `sphx_glr_auto_examples_sklearn-gridsearchcv-replacement.py` instead.
# ## Problem statement
# Tuning the hyper-parameters of a machine learning model is often carried out
# using an exhaustive exploration of (a subset of) the space all hyper-parameter
# configurations (e.g., using :obj:`sklearn.model_selection.GridSearchCV`), which
# often results in a very time consuming operation.
# In this notebook, we illustrate how to couple :class:`gp_minimize` with sklearn's
# estimators to tune hyper-parameters using sequential model-based optimisation,
# hopefully resulting in equivalent or better solutions, but within less
# evaluations.
# Note: scikit-optimize provides a dedicated interface for estimator tuning via
# :class:`BayesSearchCV` class which has a similar interface to those of
# :obj:`sklearn.model_selection.GridSearchCV`. This class uses functions of skopt to perform hyperparameter
# search efficiently. For example usage of this class, see
# `sphx_glr_auto_examples_sklearn-gridsearchcv-replacement.py`
# example notebook.

print(__doc__)
import numpy as np


# ## Objective
# To tune the hyper-parameters of our model we need to define a model,
# decide which parameters to optimize, and define the objective function
# we want to minimize.

from sklearn.datasets import load_boston
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score

boston = load_boston()
X, y = boston.data, boston.target
n_features = X.shape[1]

# gradient boosted trees tend to do well on problems like this
reg = GradientBoostingRegressor(n_estimators=50, random_state=0)


# Next, we need to define the bounds of the dimensions of the search space
# we want to explore and pick the objective. In this case the cross-validation
# mean absolute error of a gradient boosting regressor over the Boston
# dataset, as a function of its hyper-parameters.


from skopt.space import Real, Integer
from skopt.utils import use_named_args


# The list of hyper-parameters we want to optimize. For each one we define the
# bounds, the corresponding scikit-learn parameter name, as well as how to
# sample values from that dimension (`'log-uniform'` for the learning rate)
space  = [Integer(1, 5, name='max_depth'),
          Real(10**-5, 10**0, "log-uniform", name='learning_rate'),
          Integer(1, n_features, name='max_features'),
          Integer(2, 100, name='min_samples_split'),
          Integer(1, 100, name='min_samples_leaf')]

# this decorator allows your objective function to receive a the parameters as
# keyword arguments. This is particularly convenient when you want to set
# scikit-learn estimator parameters
@use_named_args(space)
def objective(**params):
    reg.set_params(**params)

    return -np.mean(cross_val_score(reg, X, y, cv=5, n_jobs=-1,
                                    scoring="neg_mean_absolute_error"))


# ## Optimize all the things!
# With these two pieces, we are now ready for sequential model-based
# optimisation. Here we use gaussian process-based optimisation.

from skopt import gp_minimize
res_gp = gp_minimize(objective, space, n_calls=50, random_state=0)

"Best score=%.4f" % res_gp.fun

# In[11]:

print("""Best parameters:
- max_depth=%d
- learning_rate=%.6f
- max_features=%d
- min_samples_split=%d
- min_samples_leaf=%d""" % (res_gp.x[0], res_gp.x[1],
                            res_gp.x[2], res_gp.x[3],
                            res_gp.x[4]))


# ## Convergence plot

from skopt.plots import plot_convergence
plot_convergence(res_gp)

-------------------------------------------------------------------------------------------------------------------

Best parameters:
- max_depth=5
- learning_rate=0.123703
- max_features=13
- min_samples_split=100
- min_samples_leaf=1
----------------------------------------------------------------------------------------------------
베이지안 옵티마이제이션을 CNN 딥러닝에 활용한 예이다.
딥러닝의 초월매개변수들을 최적화한다.

# architecture of CNN
def buildModel(neurons,drop,nb_batch,hidden_layers):
    nb_epoch=100
    neurons=int(neurons)
    hidden_layers=int(hidden_layers)
#     input layer
    inputs = Input(shape=x_train.shape[1:])
    x = Conv2D(neurons, (3, 3), padding='same', activation='relu')(inputs)
    x = MaxPooling2D((2, 2), padding='same')(x)
#     hidden layers
    if hidden_layers !=0:
        for i in range(1,hidden_layers+1):
            x = Conv2D(neurons*(2**i), (3, 3), padding='same', activation='relu')(x)
            x = MaxPooling2D((2, 2), padding='same')(x)
    x = Flatten()(x)
    x = Dense(neurons*(2**(hidden_layers+1)), activation='relu')(x)
    x = Dropout(drop)(x)
#    output
    predictions = Dense(nb_classes, activation='softmax')(x)
#     modeling
    model = Model(inputs=inputs, outputs=predictions)
    model.compile(optimizer='rmsprop', loss='binary_crossentropy', metrics=['accuracy'])
    early_stopping = EarlyStopping(patience=10, verbose=1)
#     learning
    hist = model.fit(x_train, y_train, validation_split=0.2, epochs=int(nb_epoch), batch_size=int(nb_batch),callbacks=[early_stopping])
#     evaluate the model with test data
    score = model.evaluate(x_test, y_test)
    print("loss : ",score[0])
    print("accuracy : ",score[1])
    acc.append(score[1])
    return score[0]*(-1)

#     baeysian optimization
def bayesOpt():
    pbounds = {
#     define the range of each parameter
    'drop' : (0,0.9),
"    neurons" : (2,50),
    'nb_batch' : (4,100),
    'hidden_layers' : (0,6)
    }

    optimizer = BayesianOptimization(f=buildModel, pbounds=pbounds)

#    init_points: number of initial points
#    n_iter: number of iteration
#    acq:acquire function, ucb as default
    optimizer.maximize(init_points=10, n_iter=30, acq='ucb')
    return optimizer

#     run bayesian optimization
acc =[]
result = bayesOpt()
result.res
result.max

#     save results
all_params=[]
target=[]
for i in range(len(result.res)):
    all_params.append(result.res[i]["params"])
    target.append(result.res[i]["target"])

print("all_params:")
print(all_params)
print("target:")
print(target)

params_name=list(all_params[0].keys())

with open(path+'/'+title+'.csv', 'w') as f:
    label=["acc","target(loss)"]+params_name
    writer = csv.writer(f)
    writer.writerow(label)
    for i in range(len(all_params)):      
        eachparams=[acc[i],target[i]]+list(all_params[i].values())
        writer.writerow(eachparams)
---------------------------------------------------------------------------------------------------------------------
아래의 프로그램을 보면 기존의 데이터와 새로 생성된 데이터가 기존의 데이테에 어떻게 순차적으로 추가되는지를 알 수 있다.
물론, 새로운 데이터 값을 얻어내는 것은 매우 값비싼 계산이라는 점을 염두에 두어야 한다.
값싼 계산으로서 대리모델을 이용한다. 대리모델은 하나의 함수이다. 대리 모델은 근사 함수이다. 
기계학습으로 함수를 만들어낸 것이다.
현재 상황에서 대리모델을 이용해서, 근사적으로 최고의 값을 줄 것으로 예측되는 위치를 찾아낸다.
찾아낸, 예측한 위치에서 소위 값비싼 계산을 수행한다.  
이 데이터는 기존의 데이터에 추가될 것이다.
데이터가 갱신되었기 때문에, 대리모델은 다시 fit을 해야할 것이다. 


# example of bayesian optimization for a 1d function from scratch
from math import sin
from math import pi
from numpy import arange
from numpy import vstack
from numpy import argmax
from numpy import asarray
from numpy.random import normal
from numpy.random import random
from scipy.stats import norm
from sklearn.gaussian_process import GaussianProcessRegressor
from warnings import catch_warnings
from warnings import simplefilter
from matplotlib import pyplot

# objective function
def objective(x, noise=0.1):
    noise = normal(loc=0, scale=noise)
    return (x**2 * sin(5 * pi * x)**6.0) + noise

# surrogate or approximation for the objective function
def surrogate(model, X):
    # catch any warning generated when making a prediction
    with catch_warnings():
        # ignore generated warnings
        simplefilter("ignore")
        return model.predict(X, return_std=True)

# probability of improvement acquisition function
def acquisition(X, Xsamples, model):
    # calculate the best surrogate score found so far
    yhat, _ = surrogate(model, X)
    best = max(yhat)
    # calculate mean and stdev via surrogate function
    mu, std = surrogate(model, Xsamples)
    mu = mu[:, 0]
    # calculate the probability of improvement
    probs = norm.cdf((mu - best) / (std+1E-9))
    return probs

# optimize the acquisition function
def opt_acquisition(X, y, model):
    # random search, generate random samples
    Xsamples = random(100)
    Xsamples = Xsamples.reshape(len(Xsamples), 1)
    # calculate the acquisition function for each sample
    scores = acquisition(X, Xsamples, model)
    # locate the index of the largest scores
    ix = argmax(scores)
    return Xsamples[ix, 0]

# plot real observations vs surrogate function
def plot(X, y, model):
    # scatter plot of inputs and real objective function
    pyplot.scatter(X, y)
    # line plot of surrogate function across domain
    Xsamples = asarray(arange(0, 1, 0.001))
    Xsamples = Xsamples.reshape(len(Xsamples), 1)
    ysamples, _ = surrogate(model, Xsamples)
    pyplot.plot(Xsamples, ysamples)
    # show the plot
    pyplot.show()

# sample the domain sparsely with noise
X = random(100)
y = asarray([objective(x) for x in X])
# reshape into rows and cols
X = X.reshape(len(X), 1)
y = y.reshape(len(y), 1)
# define the model
model = GaussianProcessRegressor()
# fit the model
model.fit(X, y)
# plot before hand
plot(X, y, model)
# perform the optimization process
for i in range(100):
    # select the next point to sample
    x = opt_acquisition(X, y, model)
    # sample the point
    actual = objective(x)
    # summarize the finding
    est, _ = surrogate(model, [[x]])
    print('>x=%.3f, f()=%3f, actual=%.3f' % (x, est, actual))
    # add the data to the dataset
    X = vstack((X, [[x]]))
    y = vstack((y, [[actual]]))
    # update the model
    model.fit(X, y)

# plot all samples and the final surrogate function
plot(X, y, model)
# best result
ix = argmax(y)
print('Best Result: x=%.3f, y=%.3f' % (X[ix], y[ix]))





핑백

덧글

  • 바죠 2020/02/11 17:13 # 답글

    def black_box_function(x, y):
    """Function with unknown internals we wish to maximize.

    This is just serving as an example, for all intents and
    purposes think of the internals of this function, i.e.: the process
    which generates its output values, as unknown.
    """
    return -x ** 2 - (y - 1) ** 2 + 1

    from bayes_opt import BayesianOptimization

    # Bounded region of parameter space
    pbounds = {'x': (2, 4), 'y': (-3, 3)}

    optimizer = BayesianOptimization(
    f=black_box_function,
    pbounds=pbounds,
    random_state=1,
    )



  • 바죠 2020/02/11 17:18 # 답글

    https://machinelearningmastery.com/what-is-bayesian-optimization/
  • 바죠 2020/02/12 10:52 # 답글

    http://krasserm.github.io/2018/03/19/gaussian-processes/
    https://github.com/krasserm/bayesian-machine-learning/blob/master/gaussian_processes.ipynb
    http://krasserm.github.io/2018/03/21/bayesian-optimization/
  • 바죠 2021/02/28 13:52 # 답글

    def fun(x1,x2):
    return -((x1-1)**2+ (x2-2.5)**2)
    #BO = BayesianOptimization(fun, {'x1': (0., 2.), 'x2': (0., 3.) }, verbose=2)
    list1=['x1','x2']
    list2=[(0. , 2.), (0., 3.)]
    dict( zip( list1, list2))
    BO = BayesianOptimization(fun, dict( zip( list1, list2)) , verbose=2)
    BO.maximize(init_points=2, n_iter=20)
    for i, res in enumerate(BO.res):
    print("Iteration {}: nt{}".format(i, res))
    print(BO.max)
  • 바죠 2021/07/30 11:44 # 답글

    https://news.naver.com/main/read.naver?mode=LSD&mid=shm&sid1=105&oid=032&aid=0003088775
댓글 입력 영역

최근 포토로그



MathJax