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

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

Intelligence is to find the most probable hypothesis for given training data.


H: 가설
D: 데이터

베이즈 추론(Bayesian inference)은 통계적 추론의 한 방법이다.
추론 대상의 사전 확률과 추가적인 정보를 통해 해당 대상의 사후 확률을 추론하는 방법이다. 
베이즈 추론은 베이즈 확률론을 기반으로 하며, 이는 추론하는 대상을 확률변수로 보아 그 변수의 확률분포를 추정하는 것을 의미한다.

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

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

베이즈 확률론 자체가 기계학습을 내포하고 있다. 추가되는 정보를 순차적으로 고려하면 체계적으로 본질적 실체에 대한 해석을 해 낼 수 있다. 단번에 모든것이 정리되고 성립되는 것이 아니다. 
모든 결과물에슨 추가적인 노력이 있어야 보다 더 정확한 해석이 가능하다.
지금까지 정보를 모두 종합하고 앞으로의 실험을 준비하는 것이다.


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

cos(x)라는 함수를 생각해 보자. 
여기서 우리는 cos(x)를 예로 들었다. 
우리의 가정하나를 잊지 말아야 한다. 이 함수 계산이 엄청나게 시간이 많이 걸린다고 가정하고 있다.
유한한 개수의 실측을 통해서 얻을 수 있는 함수값들을 생각한다. 
당연히 함수꼴에 대한 한정적인 정보만을 취할 수밖에 없다.
적절히 분포되어 있는 실측 함수값들이 있다고 가정한다.
4, 8, 16 points.
이러한 조건들에 대한 확률 분포를 얻어낸다. 조건부 확률을 활용한다. 평균, 편차만을 추구한다. 정규분포를 가정한다. 
실제 측정된 데이터를 기반으로 한다. 실측된 데이터가 많을수록 정확할 것이다.
확률분포에서 평균과 편차를 얻어낸다.
나아가서, 기대할 수 있는 새로운 점을 찾아낸다. (예측한다.) 추천한다.
추전받은 점에서 실측을 수행한다.
실측된 데이터는 데이터로 추가되게 된다.
주어진 데이터로부터, 이들에 합당한 조건부 확률을 표현할 수 있는 새로운 평균, 편차를 찾아낸다.
평균, 편차만 알고 있으면 추론을 할 수 있다고 가정한다.
이것으로 부터 새로운 추론이 가능하다. 새로운 가능성이 있는 영역을 간접적으로 확인할 수 있다.



기계학습에서 사용되는 모델의 성능 측정에서 cross validation score 함수는 연속함수일리가 없다. 이 함수 계산은 제법 시간이 걸리는 계산이다. 데이터가 많을 수록 더욱더 오랜 시간 계산을 수행할 것이다.
인공신경망에서도 몇 가지 변수들을 초월매개변수로 생각할 수 있다. 또한, 모델을 함수로 만들 수 있다. 즉, 바꾸고자 하는 몇 가지를 변수로 만들어서 모델을 정의할 수 있다. 예를 for loop를 이용하여 hidden layer수를 변수로 만들 수 있다. 
이 변수 또한 Bayesian optimization으로 최적화 할 수 있다. CNN의 경우도 마찬가지이다. 몇 가지를 변수로 만드는 것이 가능하다.
유한한 데이터 기반의 정밀도 측정은 불연속 함수를 사실상 가정하고 있다.
함수의 미분은 해석적으로 얻기 어렵다. 당연하다. 

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



베이지안 옵티마이제이션(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을 해야할 것이다. 
결코 병렬화 될 수 없는 알고리듬이다.
---------------------------------------------------------------------------------------------------------------------


핑백

덧글

  • 바죠 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)
  • 바죠 2022/06/19 09:46 # 답글

    # 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]))
  • 바죠 2022/06/19 09:45 # 답글

    https://machinelearningmastery.com/what-is-bayesian-optimization/
댓글 입력 영역

최근 포토로그