'Returns log-marginal likelihood of theta for training data.
Parameters
theta : array-like, shape = (n_kernel_params,) or None
Kernel hyperparameters for which the log-marginal likelihood is
evaluated. If None, the precomputed log_marginal_likelihood
of ``self.kernel_.theta`` is returned.
eval_gradient : bool, default: False
If True, the gradient of the log-marginal likelihood with respect
to the kernel hyperparameters at position theta is returned
additionally. If True, theta must not be None.
Returns
log_likelihood : float
Log-marginal likelihood of theta for training data.
log_likelihood_gradient : array, shape = (n_kernel_params,), optional
Gradient of the log-marginal likelihood with respect to the kernel
hyperparameters at position theta.
Only returned when eval_gradient is True.'

def log_marginal_likelihood(self, theta=None, eval_gradient=False):

    if (theta is None):
        if eval_gradient:
            raise ValueError('Gradient    can    only    be    evaluated    for    theta!=None')
        return self.log_marginal_likelihood_value_
    kernel = self.kernel_.clone_with_theta(theta)
    if eval_gradient:
        (K, K_gradient) = kernel(self.X_train_, eval_gradient=True)
    else:
        K = kernel(self.X_train_)
    K[np.diag_indices_from(K)] += self.alpha
    try:
        L = cholesky(K, lower=True)
    except np.linalg.LinAlgError:
        return (((- np.inf), np.zeros_like(theta)) if eval_gradient else (- np.inf))
    y_train = self.y_train_
    if (y_train.ndim == 1):
        y_train = y_train[:, np.newaxis]
    alpha = cho_solve((L, True), y_train)
    log_likelihood_dims = ((-0.5) * np.einsum('ik,ik->k', y_train, alpha))
    log_likelihood_dims -= np.log(np.diag(L)).sum()
    log_likelihood_dims -= ((K.shape[0] / 2) * np.log((2 * np.pi)))
    log_likelihood = log_likelihood_dims.sum((-1))
    if eval_gradient:
        tmp = np.einsum('ik,jk->ijk', alpha, alpha)
        tmp -= cho_solve((L, True), np.eye(K.shape[0]))[:, :, np.newaxis]
        log_likelihood_gradient_dims = (0.5 * np.einsum('ijl,ijk->kl', tmp, K_gradient))
        log_likelihood_gradient = log_likelihood_gradient_dims.sum((-1))
    if eval_gradient:
        return (log_likelihood, log_likelihood_gradient)
    else:
        return log_likelihood
