Maybe it doesn’t Matern

Calling all bored and crazy Stan coders. I need a long term coding partner.

I’ve been doing a lot of completely solo work for a couple years now, and you can guess the result. I have a couple of really awesome models that I am proud of, and yet they’re (1) not quite what I want, (2) probably way more complicated and messy than necessary, and (3) sitting there, waiting for me to figure out how to be a real modeler and justify them.

One of the things I was recently working on and am a bit stuck on is a Stan implementation of the circular Matern covariance function as described in “Covariance functions for mean square differentiable processes on spheres” (Guinness & Fuentes 2013). I just learned that the quadratic exponential covariance that I had been using is technically not valid for spherical distances, and that the linear exponential covariance is non-differentiable (bad for parameters in Stan, yes?). So this would be a nice addition. And I’m pretty sure I translated the math correctly, but it doesn’t… work. A few iterations of ADVI and it diverges when using my circular Matern implementation, whereas the version of my model with the technically invalid quadratic exponential runs just fine. So maybe I’m trying to solve a problem that doesn’t exist – but I want to make it work anyway!

Here’s the code I’ve got. I’ve filtered out a bunch of unnecessary details but haven’t cleaned the remainder up…

functions {
    real ff(int k, int j) {
    vector circular_matern(vector d, int n, real alpha, matrix ffKJ, matrix chooseRJ) {
        real ap = alpha * pi();
        row_vector[2] csap = [cosh(ap), sinh(ap)];
        matrix[n-1,n] H;
        real annm1 = inv((-2 * square(alpha))^(n-1) * tgamma(n));
        vector[n] a;
        vector[rows(d)] adp = alpha * (d-pi());
        matrix[rows(d),2] csadp = append_col(cosh(adp), sinh(adp));
        vector[rows(d)] cov = zeros_vector(rows(d));
        for(k in 0:(n-1)) {
          for(r in 0:(n-2)) {
            H[r+1,k+1] = 0;
            for(j in 0:(2*r+1)) {
              if(j <= k) {
                 += chooseRJ[r+1,j+1]
                    * ffKJ[k+1,j+1]
                    * ap^(k-j)
                    * csap[((k-j+1) % 2) + 1];
        a = append_row(-annm1 * (H[,1:(n-1)] \ H[,n-1]), annm1);
        for(k in 0:(n-1)) {
          cov += a[k+1] * adp^k .* csadp[,(k % 2) + 1];
        return(cov / (2*alpha*csap[2]));
    matrix fill_sym(vector lt, int N, real c) {
        matrix[N,N] s_mat;
        int iter = 1;
        for(j in 1:(N-1)) {
            s_mat[j,j] = c;
            for(i in (j+1):N) {
                s_mat[i,j] = lt[iter];
                s_mat[j,i] = lt[iter];
                iter += 1;
        s_mat[N,N] = c;
data {
    int M[DRC];
    vector[choose(M[size(M)],2)] distSites; // distm(latlon)[lower.tri(distm(latlon)]
    int ms; // Matern smoothness parameter
transformed data {
    matrix[ms-1, 2 * ms] chooseRJ;
    matrix[ms, 2 * ms] ffKJ;
parameters {
    matrix<lower=0>[DRC+D,K] weight_scales;
    vector<lower=0>[K] siteDecay; //alpha
    vector<lower=0, upper=1>[K] siteProp; //nugget
transformed parameters {
    for(k in 1:K) {
            = fill_sym(siteProp[k]
                       * square(weight_scales[DRC,k])
                       * circular_matern(distSites, ms, inv(siteDecay[k]), ffKJ, chooseRJ),
                       square(weight_scales[DRC,k]) + 1e-10);
        // below is the alternative, which works
        //    = fill_sym(siteProp[k]
        //               * square(weight_scales[DRC,k])
        //               * exp(-square(distSites) / (2 * square(siteDecay[k]))),
        //               M[size(M)],
        //               square(weight_scales[DRC,k]) + 1e-10);
model {
    weight_scales ~ student_t(2,0,2.5);
    siteDecay ~ inv_gamma(5,5);

Mer de Microbes

When I started my first blog, I initially used this title for it. I had recently begun working in French Polynesia, was helping a French postdoc, and was interested in learning the language, and I thought, “‘Sea of Microbes’ in French, that’s nice.” I created the page and started working on the blog. And then that French postdoc saw the url: ‘merdemicrobes’.

‘What?? You can’t really call it that…’

His first parsing of the spaceless text was ‘Merde Microbes’. Microbe shit.

This was hilarious to me, but I was using a university-hosted blogging platform and ultimately wanted (at least some of) my writing to be accessible to kids, so, chagrined, I came up with a more bland name and commenced my writing and sciencing adventure under that.

Since then, my adventures have seen some wild ups and downs. My blog expanded for a bit with some features and posts that I was proud of. And then I got busy with other things, discovered that most of the features I was developing had already been done (better) by others, lost my motivation, and it trailed off. This blog arc largely paralleled my sentiments toward my actual science, as well. The more you learn, the more you learn that others have already learned it.

I’ve had a few years of the real world of academia and politics to disabuse myself of perfectionism and idealism. I’ve realized that I’ve tried to hold myself to standards of communication that are impossible for me to achieve, and that my constant self-censorship has led to extreme frustration, burnout, and anxiety. Though I was trying to be the best scientist and communicator possible, I ultimately became just completely stuck. So today I’m trying something new. I’m laying claim to ‘Mer de Microbes’ again. Fully embracing its original intent, but also the accidental double-meaning, the grammatical imperfection, the fact that I’m not French, and everything else that will come with the name. It may change again later. But for now, I am going to try and recommence my journey through a sea of microbes. Or through microbe shit. It doesn’t really matter.