function [logL, logL_theta] = sqexp_marginal_likelihood(x, y, par)
% MARGINAL_LIKELIHOOD Calculate the marginal likelihood (eqn (5.8) of
% Rasmussen and Williams) and (if necessary) its derivative (eqn (5.9)).
% This can be used to optimise the hyper-parameters of the GP.

% Calculate the covariance matrix with derivatives
[K, K_sigma_n, K_sigma_f, K_l] = sqexp_covar(x, x, par);

% Calculate the Cholesky decomposition of the covariance matrix
[n, ~] = size(K);
R = chol(K); % K = R'*R

% Calculate the log determinant of the covariance matrix
logdet_K = sum(2*log(R(1:n+1:end)));

% Calculate the log likelihood
alpha = R \ (R' \ y');
logL = -0.5*y*alpha - 0.5*logdet_K - 0.5*n*log(2*pi);

% Calculate the derivatives
if nargout > 1
    inner = (alpha*alpha' - R \ (R' \ eye(n)))';
    logL_theta = zeros(size(par));
    logL_theta(1) = 0.5*sum(sum(inner.*K_sigma_n));
    logL_theta(2) = 0.5*sum(sum(inner.*K_sigma_f));
    l = par(3:end);
    if length(l) == 1
        logL_theta(3) = 0.5*sum(sum(inner.*K_l, 2));
    else
        for i = 1:length(l)
            logL_theta(i + 2) = 0.5*sum(sum(inner.*K_l(:, :, i)));
        end
    end
end
