가장 쉬운 가우시인 프로세스, 기대 향상, 그리고 베이지안 옵티마이제이션 by 바죠

가장 쉬운 가우시인 프로세스, 기대 향상, 그리고 베이지안 옵티마이제이션 (Bayesian optimization)


가우시안 프로세스 (regression): 회귀
가우시안 프로세스는 베이지안 접근 방식을 사용하여 데이터에 적합한 모든 가능한 함수 f (x)에 대한 확률 분포를 추론한다.
모델은 비모수 커널 기반의 확률적 모델이다.
모델이 확률적이므로 훈련된 모델을 사용하여 예측 구간을 계산한다.

어떻게 사전 정보를 통합적으로 관리하는가? 
이는 공분산 행렬을 결정하는 커널을 통해서 사전의 정보를 통합적으로 관리한다.

관측들 사이의 유사성은 가우스 분포로 모델링되며 인접 이웃과의 공동 확률 분포를 갖는다.
이러한 의미에서 각 관측치는 다차원에서 좌표를 갖는다.
추가 차원을 포함하도록 지수 커널을 확장 할 수도 있다.
커널 또는 공분산 함수는 모든 x 및 x ' 관측치 쌍 사이의 유사성을 정의한다. 
커널은 일반적인 발전한 내적으로 볼 수 있다.

새로운 위치에서의 예측은 지나간 관측에서 얻어낸 데이터의 선형 조합니다.
평균 예측(A) 그리고 분산 예측(B)을 알면 새로운 위치에서의 추정을 실행할 수 있다.
새로운 위치
\[x^* \]

과거 데이터
\[ \{x\} \]
\[ \{f(x)\} \]

새 위치와 과거 위치 상관관계 :
\[k(\theta, x^*, x) \]
\[ p(y^* \mid x^*, f(x), x)  = N(y^*\mid A, B) \] 
\[ m(x)=E[f(x)] \]
\[ k(x,x')=E[(f(x)-m(x))(f(x')-m(x'))]  \]
\[ f(x^*) \sim GP(m(x),k(x,x')) \]

입력 벡터 x에 대해서 특징 벡터를 활용한다. 
\[ \phi(x) \]
\[ f(x) = \phi^T(x) N(0, \Sigma_p)\]
\[ k(x,x') = \phi(x)\Sigma_p \phi(x') \]
\[ f^* \mid f \sim N(...) \]

이렇게 얻어진 데이터를 이용하여 회귀함수를 만든다.
이것이 가우시안 프로세스 방법이다. 


베이지안 옵티마이제이션 (Bayesian optimization): 최적화 방법  http://incredible.egloos.com/7479039

베이즈 확률론(Bayesian probability)은 확률을 '지식 또는 믿음의 정도를 나타내는 양'으로 해석하는 확률론이다.
확률을 발생 빈도(frequency)나 어떤 시스템의 물리적 속성으로 여기는 것과는 다른 해석이다. 이 분야의 선구자인 18세기 통계학자 토머스 베이즈의 이름을 따서 명명되었다.
어떤 가설의 확률을 평가하기 위해서 사전 확률을 먼저 밝히고 새로운 관련 데이터에 의한 새로운 확률값을 변경한다.
데이터가 주어지기 전에 이미 어느 정도 확률값을 예측하고 있을 때 이를 새로 수집한 데이터와 합쳐서 최종 결과에 반영할 수 있다.

인공지능 분야에서 잘 활용할 수 있는 통계학이다.
데이터가 많을 수록 더 정교한 예측이 가능하다.


가우시안 프로세스로부터 얻어낸 회귀함수를 활용한다. 
좋은 목적함수를 줄 것으로 예상되는 점을 추정한다. 오차까지 함께 추정한다.
지금까지 얻어낸 것 보다 더 좋은 해가 있는 곳을 예측한다.
회귀함수를 활용한 예측을 수행한다. 새로운 점에서의 목적함수 값을 예측한다.
가볍게 계산할 수 있고 간단한 최적화를 할 수 있다.

기대 향상 (expected improvement) : 현재 최고값보다 더 나아진 정도, 기대되는 개선의 정도 
f(x+) :  현재까지 최고값.
x+   :  현재까지 최고의 해
xi = 0.01 로 고정한다. 대리 모델을 활용한 기대 향상의 예측을 생각해 볼 수 있다. 

\[ f(\vec{x^{+}}) \]
\[ \vec{x^{+}} \]
\[ \xi=0.01 \]

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

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

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

제안된 점에서 목적함수 계산을 수행한다.
이 정보를 기존의 데이터에 추가시킨다. 이제 데이터의 양이 증가한 상황이 된다.
처음과 마찬가지로 보강된 데이터를 이용하여 새로운 점을 예측하는 방식으로 계산을 반복한다.

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





여기에서 사용하는 함수 두 개가 있다. 모두 정규분포와 관련된 함수들이다. 

cdf : 주어진 확률 변수가 특정 값보다 작거나 같은 확률을 나타내는 함수이다.
Φ and ϕ are the CDF and PDF of the standard normal distribution, respectively.
첫항: exploitation
둘째항: exploration
higher /www.w3.org/1998/Math/MathML"">ξ" role="presentation" style="display: inline; line-height: normal; font-size: medium; overflow-wrap: normal; white-space: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; padding: 0px; margin: 0px; font-family: "Malgun Gothic"; position: relative;">ξξ values lead to more exploration

#   http://krasserm.github.io/2018/03/21/bayesian-optimization/
#   Modified by In-Ho Lee, KRISS, February 14, 2020.
import numpy as np
from scipy.stats import norm
def expected_improvement(X,X_sample,Y_sample,gpr,xi=0.01):
    '''
    Computes the EI at points X based on existing samples X_sample
    and Y_sample using a Gaussian process surrogate model.
    Args:
        X: Points at which EI shall be computed (m x d).
        X_sample: Sample locations (n x d).
        Y_sample: Sample values (n x 1).
        gpr: A GaussianProcessRegressor fitted to samples.
        xi: Exploitation-exploration trade-off parameter.
    Returns:
        Expected improvements at points X.
    '''
    mu, sigma=gpr.predict(X,return_std=True)
    sigma1=np.zeros_like(sigma)
    sigma1[:]=sigma1[:]+1e-16
    sigma=sigma+sigma1
    mu_sample=gpr.predict(X_sample)
    sigma = sigma.reshape(-1, 1)
#   Needed for noise-based model, otherwise use np.max(Y_sample). See also section 2.4 in [...]
    mu_sample_opt = np.max(mu_sample)
    with np.errstate(divide='warn'):
        imp = mu - mu_sample_opt - xi
        Z = imp / sigma
        ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
        ei[sigma == 0.0] = 0.0
    return ei
from scipy.optimize import minimize
def propose_location(acquisition,X_sample,Y_sample,gpr,bounds, n_restarts=300):
    '''
    Proposes the next sampling point by optimizing the acquisition function.
    Args:
        acquisition: Acquisition function.
        X_sample: Sample locations (n x d).
        Y_sample: Sample values (n x 1).
        gpr: A GaussianProcessRegressor fitted to samples.
    Returns:
        Location of the acquisition function maximum.
    '''
    ndim = X_sample.shape[1]
    min_val = 1
    min_x = None
    def min_obj(X):
#   Minimization objective is the negative acquisition function
        return -acquisition(X.reshape(-1,ndim),X_sample,Y_sample,gpr)
#   Find the best optimum by starting from n_restart different random points.
    for x0 in np.random.uniform(bounds[:, 0],bounds[:, 1], size=(n_restarts, ndim)):
        res = minimize(min_obj,x0=x0,bounds=bounds,method='L-BFGS-B')
        if res.fun < min_val:
            min_val = res.fun
            min_x = res.x
    return min_x.reshape(-1,ndim)


ndim=5
ninit=10*ndim
bounds=np.zeros((ndim,2))
for i in range(ndim):
    for j in range(2):
        if j == 0 :
           bounds[i,j]= -1.00+float(i)
           bounds[i,j]= -3.00
        if j == 1 :
           bounds[i,j]= 1.00+float(i)
           bounds[i,j]= 3.00
noise = 1e-8
def f(X,noise=noise):
    y=np.zeros((X.shape[0],))
    ndim = X.shape[1]
    if False:
        for i in range(X.shape[0]):
            y[i]=-10.*ndim
            for j in range(X.shape[1]):
                y[i]=y[i]-X[i,j]**2+10.*np.cos(2.*(np.pi)*X[i,j])
    if True:
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                y[i]=y[i]+np.sin(X[i,j])
    if False:
        for i in range(X.shape[0]):
            for j in range(X.shape[1]):
                y[i]=y[i]-(X[i,j]-float(j))**2
    return y
X_init=np.zeros((ninit,ndim))
for i in range(ninit):
    for j in range(ndim):
        X_init[i,j]=bounds[j,0]+(bounds[j,1]-bounds[j,0])*np.random.random()
Y_init=f(X_init)
for i in range(ninit):
    print(X_init[i,:],Y_init[i])
best_obj=Y_init[0]
best_sol=X_init[0,:]
for i in range(ninit):
    if best_obj < Y_init[i]:
       best_obj=Y_init[i]
       best_sol=X_init[i,:]
       print(best_obj)
from sklearn.gaussian_process import  GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernelMatern
#   Gaussian process with Matern kernel as surrogate model
m52=ConstantKernel(1.0) * Matern(length_scale=1.0, nu=2.5)
gpr=GaussianProcessRegressor(kernel=m52, alpha=noise**2)
X_sample=X_init
Y_sample=Y_init
n_iter =50*ndim
for i in range(n_iter):
    gpr.fit(X_sample,Y_sample)
#   Obtain next sampling point from the acquisition function (expected_improvement)
    X_next=propose_location(expected_improvement,X_sample,Y_sample,gpr,bounds,n_restarts=30*ndim)
    Y_next=f(X_next,noise)
#   print(X_next.ravel(),Y_next)
    if best_obj < Y_next:
       best_obj=Y_next
       best_sol=X_next
       print(X_next.ravel(),Y_next)
    X_sample=np.vstack((X_sample,X_next))
    Y_sample=np.hstack((Y_sample,Y_next))
print(best_sol.ravel(),best_obj)
import matplotlib.pyplot as plt
plt.rc('font', size=14)          # controls default text sizes
#plt.rc('figure', titlesize=14)  # fontsize of the figure title
#plt.rc('axes', titlesize=14)     # fontsize of the axes title
#plt.rc('axes', labelsize=18)    # fontsize of the x and y labels
#plt.rc('xtick', labelsize=16)    # fontsize of the tick labels
#plt.rc('ytick', labelsize=16)    # fontsize of the tick labels
#plt.rc('legend', fontsize=14)    # legend fontsize
x=np.zeros((Y_sample.shape[0],))
for i in range(Y_sample.shape[0]):
    x[i]=i+1
fig=plt.figure()
plt.plot(x[:ninit], Y_sample[:ninit], marker='+', color='red', linewidth=2, markersize=8, label='Random')
plt.plot(x[ninit:], Y_sample[ninit:], marker='o', color='blue', linewidth=2, markersize=8, label='Gaussian Process')
plt.xlabel('Samples', fontsize=20)
plt.ylabel('Objective function', fontsize=20)
plt.legend()
fig.tight_layout()
fig.savefig('gpsamples.eps')
#plt.show()


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

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


#!/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

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


Imagination is more important than knowledge (지식보다 상상) : 가장 쉬운 베이지안 옵티마이제이션 [Python] (egloos.com)

---------------------------------------------------------------------------------------------------------------------
베이즈 추정에는 축차합리성이 있다.

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


# 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]))

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


핑백

덧글

  • 2021/01/14 09:11 # 답글 비공개

    비공개 덧글입니다.
댓글 입력 영역

최근 포토로그



MathJax