<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Bayesian optimization
This notebook was made with the following version of george
Step1:  In this tutorial, we'll show a very simple example of implementing "Bayesian optimization" using george.
Now's not the time to get into a discussion of the issues with the name given to these methods, but I think that the "Bayesian" part of the title comes from the fact that the method relies on the (prior) assumption that the objective function is smooth.
The basic idea is that you can reduce the number of function evaluations needed to minimize a black-box function by using a GP as a surrogate model.
This can be huge if evaluating your model is computationally expensive.
I think that the classic reference is Jones et al. (1998) and the example here will look a bit like their section 4.1.
First, we'll define the scalar objective that we want to minimize in the range $-5 \le x \le 5$.
Step2:  Now, for the "Bayesian" optimization, the basic procedure that we'll follow is
<ASSISTANT_TASK:>
Python Code:
import george
george.__version__
Explanation: Bayesian optimization
This notebook was made with the following version of george:
End of explanation
import numpy as np
import matplotlib.pyplot as plt
def objective(theta):
    return -0.5 * np.exp(-0.5*(theta - 2)**2) - 0.5 * np.exp(-0.5 * (theta + 2.1)**2 / 5) + 0.3
t = np.linspace(-5, 5, 5000)
plt.figure(figsize=(8, 5))
plt.plot(t, objective(t))
plt.ylim(-0.37, 0.37)
plt.xlim(-5, 5)
plt.xlabel("$\\theta$")
plt.ylabel("objective");
Explanation: In this tutorial, we'll show a very simple example of implementing "Bayesian optimization" using george.
Now's not the time to get into a discussion of the issues with the name given to these methods, but I think that the "Bayesian" part of the title comes from the fact that the method relies on the (prior) assumption that the objective function is smooth.
The basic idea is that you can reduce the number of function evaluations needed to minimize a black-box function by using a GP as a surrogate model.
This can be huge if evaluating your model is computationally expensive.
I think that the classic reference is Jones et al. (1998) and the example here will look a bit like their section 4.1.
First, we'll define the scalar objective that we want to minimize in the range $-5 \le x \le 5$.
End of explanation
from george import kernels
from scipy.special import erf
from scipy.optimize import minimize
N_init = 4
train_theta = np.linspace(-5, 5, N_init + 1)[1:]
train_theta -= 0.5 * (train_theta[1] - train_theta[0])
train_f = objective(train_theta)
gp = george.GP(np.var(train_f) * kernels.Matern52Kernel(3.0),
               fit_mean=True)
gp.compute(train_theta)
def nll(params):
    gp.set_parameter_vector(params)
    g = gp.grad_log_likelihood(train_f, quiet=True)
    return -gp.log_likelihood(train_f, quiet=True), -g
fig, axes = plt.subplots(2, 2, figsize=(8, 6))
j = 0
old_min = None
converged = False
for i in range(1000):
    # Update the GP parameters
    soln = minimize(nll, gp.get_parameter_vector(), jac=True)
    # Compute the acquisition function
    mu, var = gp.predict(train_f, t, return_var=True)
    std = np.sqrt(var)    
    f_min = np.min(train_f)
    chi = (f_min - mu) / std 
    Phi = 0.5 * (1.0 + erf(chi / np.sqrt(2)))
    phi = np.exp(-0.5 * chi**2) / np.sqrt(2*np.pi*var)
    A_ei = (f_min - mu) * Phi + var * phi
    A_max = t[np.argmax(A_ei)]
    # Add a new point
    train_theta = np.append(train_theta, A_max)
    train_f = np.append(train_f, objective(train_theta[-1]))
    gp.compute(train_theta)
    
    # Estimate the minimum - I'm sure that there's a better way!
    i_min = np.argmin(mu)
    sl = slice(max(0, i_min - 1), min(len(t), i_min + 2))
    ts = t[sl]
    D = np.vander(np.arange(len(ts)).astype(float))
    w = np.linalg.solve(D, mu[sl])
    minimum = ts[0] + (ts[1] - ts[0]) * np.roots(np.polyder(w[::-1]))
    
    # Check convergence
    if i > 0 and np.abs((old_min - minimum) / minimum) < 1e-5:
        converged = True
    old_min = float(minimum[0])
    
    # Make the plots
    if converged or i in [0, 1, 2]:
        ax = axes.flat[j]
        j += 1
        ax.plot(t, objective(t))
        ax.plot(t, mu, "k")
        ax.plot(train_theta[:-1], train_f[:-1], "or")
        ax.plot(train_theta[-1], train_f[-1], "og")
        ax.fill_between(t, mu+std, mu-std, color="k", alpha=0.1)
        if i <= 3:
            ax2 = ax.twinx()
            ax2.plot(t, A_ei, "g", lw=0.75)
            ax2.set_yticks([])
        ax.axvline(old_min, color="k", lw=0.75)
        ax.set_ylim(-0.37, 0.37)
        ax.set_xlim(-5, 5)
        ax.set_yticklabels([])
        ax.annotate("step {0}; {1:.3f}".format(i+1, old_min), xy=(0, 1),
                    xycoords="axes fraction", ha="left", va="top",
                    xytext=(5, -5), textcoords="offset points",
                    fontsize=14)
    
    if converged:
        break
plt.tight_layout()
print("{0} model evaluations".format(len(train_f)))
Explanation: Now, for the "Bayesian" optimization, the basic procedure that we'll follow is:
Start by evaluating the model at a set of points. In this case, we'll start with a uniform grid in $x$.
Fit a GP (optimize the hyperparameters) to the set of training points.
Find the input coordinate that maximizes the "expected improvement" (see Section 4 of Jones+ 1998). For simplicity, we simply use a grid search to maximize this, but this should probably be a numerical optimization in any real application of this method.
At this new coordinate, run the model and add this as a new training point.
Return to step 2 until converged. We'll judge convergence using relative changes in the location of the minimum.
End of explanation
<END_TASK>