'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_)
    (Z, (pi, W_sr, L, b, a)) = self._posterior_mode(K, return_temporaries=True)
    if (not eval_gradient):
        return Z
    d_Z = np.empty(theta.shape[0])
    R = (W_sr[:, np.newaxis] * cho_solve((L, True), np.diag(W_sr)))
    C = solve(L, (W_sr[:, np.newaxis] * K))
    s_2 = (((-0.5) * (np.diag(K) - np.einsum('ij,    ij    ->    j', C, C))) * ((pi * (1 - pi)) * (1 - (2 * pi))))
    for j in range(d_Z.shape[0]):
        C = K_gradient[:, :, j]
        s_1 = ((0.5 * a.T.dot(C).dot(a)) - (0.5 * R.T.ravel().dot(C.ravel())))
        b = C.dot((self.y_train_ - pi))
        s_3 = (b - K.dot(R.dot(b)))
        d_Z[j] = (s_1 + s_2.T.dot(s_3))
    return (Z, d_Z)
