<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A gentle introduction to Gaussian Process Regression
This notebook was made with the following version of george
Step1:  We'll start by generating some fake data (from a sinusoidal model) with error bars
Step2:  Now, we'll choose a kernel (covariance) function to model these data, assume a zero mean model, and predict the function values across the full range. The full kernel specification language is documented here but here's an example for this dataset
Step3:  The gp model provides a handler for computing the marginalized likelihood of the data under this model
Step4:  So we can use this—combined with scipy's minimize function—to fit for the maximum likelihood parameters
Step5:  And plot the maximum likelihood model
<ASSISTANT_TASK:>
Python Code:
import george
george.__version__
Explanation: A gentle introduction to Gaussian Process Regression
This notebook was made with the following version of george:
End of explanation
import numpy as np
import matplotlib.pyplot as pl
np.random.seed(1234)
x = 10 * np.sort(np.random.rand(15))
yerr = 0.2 * np.ones_like(x)
y = np.sin(x) + yerr * np.random.randn(len(x))
pl.errorbar(x, y, yerr=yerr, fmt=".k", capsize=0)
pl.xlim(0, 10)
pl.ylim(-1.45, 1.45)
pl.xlabel("x")
pl.ylabel("y");
Explanation: We'll start by generating some fake data (from a sinusoidal model) with error bars:
End of explanation
from george import kernels
kernel = np.var(y) * kernels.ExpSquaredKernel(0.5)
gp = george.GP(kernel)
gp.compute(x, yerr)
x_pred = np.linspace(0, 10, 500)
pred, pred_var = gp.predict(y, x_pred, return_var=True)
pl.fill_between(x_pred, pred - np.sqrt(pred_var), pred + np.sqrt(pred_var),
                color="k", alpha=0.2)
pl.plot(x_pred, pred, "k", lw=1.5, alpha=0.5)
pl.errorbar(x, y, yerr=yerr, fmt=".k", capsize=0)
pl.plot(x_pred, np.sin(x_pred), "--g")
pl.xlim(0, 10)
pl.ylim(-1.45, 1.45)
pl.xlabel("x")
pl.ylabel("y");
Explanation: Now, we'll choose a kernel (covariance) function to model these data, assume a zero mean model, and predict the function values across the full range. The full kernel specification language is documented here but here's an example for this dataset:
End of explanation
print("Initial ln-likelihood: {0:.2f}".format(gp.log_likelihood(y)))
Explanation: The gp model provides a handler for computing the marginalized likelihood of the data under this model:
End of explanation
from scipy.optimize import minimize
def neg_ln_like(p):
    gp.set_parameter_vector(p)
    return -gp.log_likelihood(y)
def grad_neg_ln_like(p):
    gp.set_parameter_vector(p)
    return -gp.grad_log_likelihood(y)
result = minimize(neg_ln_like, gp.get_parameter_vector(), jac=grad_neg_ln_like)
print(result)
gp.set_parameter_vector(result.x)
print("\nFinal ln-likelihood: {0:.2f}".format(gp.log_likelihood(y)))
Explanation: So we can use this—combined with scipy's minimize function—to fit for the maximum likelihood parameters:
End of explanation
pred, pred_var = gp.predict(y, x_pred, return_var=True)
pl.fill_between(x_pred, pred - np.sqrt(pred_var), pred + np.sqrt(pred_var),
                color="k", alpha=0.2)
pl.plot(x_pred, pred, "k", lw=1.5, alpha=0.5)
pl.errorbar(x, y, yerr=yerr, fmt=".k", capsize=0)
pl.plot(x_pred, np.sin(x_pred), "--g")
pl.xlim(0, 10)
pl.ylim(-1.45, 1.45)
pl.xlabel("x")
pl.ylabel("y");
Explanation: And plot the maximum likelihood model:
End of explanation
<END_TASK>