Is this the correct code for computing the laplace approximation of the marginal data density?

function [log_mdd, log_mdd_se] = laplace_mdd(prior_fun, log_likelihood, theta_mode, hessian)
% Compute the Laplace approximation of the marginal data density

% Evaluate the log prior and log likelihood at the mode
log_prior_mode = log(prior_fun(theta_mode));
log_like_mode = log_likelihood(theta_mode);

% Compute the Laplace approximation of the marginal data density
log_mdd = log_like_mode + log_prior_mode - 0.5 * log(det(-hessian));

% Compute the standard error of the log marginal data density
log_mdd_se = 
sqrt(sum(diag(-inv(hessian))));
end     

I asked chatGPT to give me code for the laplace approximation of the marginal data density and it gave me this. The log_mdd might be right but I'm very skeptical of its computation of the standard error as the sqrt of the sum of the trace of the inverse of the negative hessian.

Answer 1

[Rishav]

Hi Jacob,

The code that you provided for finding the standard error for laplace approximation is correct.

log_mdd_se = sqrt(sum(diag(-inv(hessian))))

where, inv(hessian) is the inverse of the Hessian matrix and diag() extracts the diagonal elements of the matrix.

To know whether the standard error of the log marginal data density computed using the Laplace approximation is correct or not, you can compare it to the standard error computed using other methods, such as Monte Carlo simulation or numerical integration. If the standard error of the log marginal data density computed using the Laplace approximation is close to the standard error computed using Monte Carlo simulation or numerical integration, then it is likely that the Laplace approximation is a good approximation of the posterior density.

Regards,

Rishav Saha