Making arbitrary C++ functions available to Stan


I am trying to implement a custom function in C++ that will be ultimately used to compute the likelihood (and can be used to compute its gradient) – an example of a different gaussian process kernel and function to obtain the “K” matrix is shown below.

I am trying to follow instructions on Adding a new function with known gradients, but I gather from another thread that some function templating may be required, but I’m missing quite a few details on how to actually get it working. I understand that this is more straightforwardly done in CmdStan, so I installed it for that specifically (I have used RStan and PyStan just a couple times before that in the “standard” way).

I assume it will go something like this: write the templated C++ function, include the file under some Stan-math library folder, write the Stan program using the function and compiling it with the CmdStan interface.
Bearing in mind that I am new to Stan, can handle some C++ but I’m no developer, I’d appreciate if I could get help with the steps described above. Thanks.

#include <stan/math/rev.hpp>

namespace stan {
  namespace math {

    double exp_quad_kernel(double xa, double xb, double l1, double l2) {
        double k = exp(-pow(abs(xa-xb),2)/(pow(l1,2)+pow(l2,2)));
        return k;

    MatrixXd k_matrix(VectorXd x, ArrayXd lengthscales, double sigmaf) {

    MatrixXd K(x.rows(), x.rows());
    double l1 = lengthscales[0];
    double l2 = lengthscales[1];
    for (int i=0; i<x.cols(); i++) {
        for (int j=0; j<x.cols(); j++) {
            K(i,j) = pow(sigmaf,2)*exp_quad_kernel(x(i), x(j), l1, l2);
    return K;