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;
}
}
}
```