Hi!
There is a C++ implementation of the integrated Laplace approximation, i.e. a function that returns the approximate marginal distribution (which you call “likelihood”) and its gradient, and an experimental branch to run this function in Stan (or there should be sometime this week). Here’s the prototype stan-math code.
You can pass a function that returns the prior covariance and one the log likelihood (conditional on the latent variables you want to marginalize out), and the autodiff handles the rest, including Hessian matrices and third-order derivatives. The main constraint is that any operation in the log likelihood function must support forward mode autodiff. Neither the optimization required to compute the Laplace approximation, nor the autodiff are trivial. I have a draft that derives the differentiation rules used in the Stan code. I’m guessing that by the time you finish your PhD this document will be finished and available.
I think this should get you some of the way towards implementing EM + Laplace in Stan’s C++. From reading the doc you shared, I don’t think you’ll need the gradient of the approximate marginal likelihood but if you do, it’s there.