'Fit Gaussian process regression model.
Parameters
X : array-like, shape = (n_samples, n_features)
Training data
y : array-like, shape = (n_samples, [n_output_dims])
Target values
Returns
self : returns an instance of self.'

def fit(self, X, y):

    if (self.kernel is None):
        self.kernel_ = (C(1.0, constant_value_bounds='fixed') * RBF(1.0, length_scale_bounds='fixed'))
    else:
        self.kernel_ = clone(self.kernel)
    self._rng = check_random_state(self.random_state)
    (X, y) = check_X_y(X, y, multi_output=True, y_numeric=True)
    if self.normalize_y:
        self._y_train_mean = np.mean(y, axis=0)
        y = (y - self._y_train_mean)
    else:
        self._y_train_mean = np.zeros(1)
    if (np.iterable(self.alpha) and (self.alpha.shape[0] != y.shape[0])):
        if (self.alpha.shape[0] == 1):
            self.alpha = self.alpha[0]
        else:
            raise ValueError(('alpha    must    be    a    scalar    or    an    array    with    same    number    of    entries    as    y.(%d    !=    %d)' % (self.alpha.shape[0], y.shape[0])))
    self.X_train_ = (np.copy(X) if self.copy_X_train else X)
    self.y_train_ = (np.copy(y) if self.copy_X_train else y)
    if ((self.optimizer is not None) and (self.kernel_.n_dims > 0)):
        def obj_func(theta, eval_gradient=True):
            if eval_gradient:
                (lml, grad) = self.log_marginal_likelihood(theta, eval_gradient=True)
                return ((- lml), (- grad))
            else:
                return (- self.log_marginal_likelihood(theta))
        optima = [self._constrained_optimization(obj_func, self.kernel_.theta, self.kernel_.bounds)]
        if (self.n_restarts_optimizer > 0):
            if (not np.isfinite(self.kernel_.bounds).all()):
                raise ValueError('Multiple    optimizer    restarts    (n_restarts_optimizer>0)    requires    that    all    bounds    are    finite.')
            bounds = self.kernel_.bounds
            for iteration in range(self.n_restarts_optimizer):
                theta_initial = self._rng.uniform(bounds[:, 0], bounds[:, 1])
                optima.append(self._constrained_optimization(obj_func, theta_initial, bounds))
        lml_values = list(map(itemgetter(1), optima))
        self.kernel_.theta = optima[np.argmin(lml_values)][0]
        self.log_marginal_likelihood_value_ = (- np.min(lml_values))
    else:
        self.log_marginal_likelihood_value_ = self.log_marginal_likelihood(self.kernel_.theta)
    K = self.kernel_(self.X_train_)
    K[np.diag_indices_from(K)] += self.alpha
    try:
        self.L_ = cholesky(K, lower=True)
    except np.linalg.LinAlgError as exc:
        exc.args = ((("The    kernel,    %s,    is    not    returning    a    positive    definite    matrix.    Try    gradually    increasing    the    'alpha'    parameter    of    your    GaussianProcessRegressor    estimator." % self.kernel_),) + exc.args)
        raise
    self.alpha_ = cho_solve((self.L_, True), self.y_train_)
    return self
