Implementation for bayesian optimization with PoI acquisition function for finding the minimum point of the function

BayesianSearchOptimization.py

from scipy import minimum,maximum
from scipy.optimize import minimize
import time
from sklearn.gaussian_process import GaussianProcess
from scipy.stats import norm,entropy
Bounds = [[-16,16]]
#target function
def f(x):
    y = x**2*np.sin(x)
    return y
X = np.asarray([np.arange(x[0], x[1], 1e-1) for x in Bounds]).T
y = np.asarray([f(_x) for _x in X])
plt.figure(figsize=(15,5))
plt.subplot(3,1,1)
plt.plot(X,y)
plt.title('real target values we want to search for -200')
gp = GaussianProcess()
sampling_size = 10
x_try = np.asarray([[np.random.uniform(x[0], x[1], size=1) for x in Bounds] for _ in range(sampling_size)]).T[0].T
y_try = np.asarray([f(_x) for _x in x_try])
gp.fit(x_try,y_try)
repeat_search = 10
score_random = []
score_bayes = []
def PoI_acquisition(x,gp,fmin):
    mean, variance = gp.predict(x.reshape(-1,1), eval_MSE=True)
    if variance == 0:
        print 'it wrong here'
        return 1
    else:
        z_score = (fmin - mean)/np.sqrt(variance)
        return norm.cdf(z_score)
for i in range(repeat_search):
    #random search
    _x = np.asarray([np.random.uniform(x[0], x[1]) for x in Bounds]) 
    score_random.append(f(_x))
    #bayesian search
    fmin = np.min(y_try)
    acq_min = np.Inf
    for _ in range(50):
        __x = np.asarray([np.random.uniform(x[0], x[1]) for x in Bounds])
        guess_x = minimize(lambda x: -PoI_acquisition(x,gp,fmin),[__x], bounds = Bounds, method = 'L-BFGS-B')
        if guess_x.fun < acq_min:
            print 'better guess',guess_x.fun
            best_x = guess_x
            acq_min = guess_x.fun
    print(_x,best_x.x,best_x.fun,fmin)
    _x = best_x.x
    _y = f(_x)
    score_bayes.append(_y)
    x_try = np.concatenate((x_try,[_x]),axis=0)
    y_try = np.concatenate((y_try,[_y]),axis=0)
    gp.fit(x_try,y_try)
plt.subplot(3,1,2)
plt.plot(score_random)
plt.title('optimize by random search')
plt.subplot(3,1,3)
plt.plot(score_bayes)
plt.title('optimize by bayes search')
#okay now you can see that bayes search return better result than random search