Puzzling error when fitting multivariate normal model

Modeling
1

Hello,

I ran into a puzzling error when testing a multivariate normal model and was hoping you could help me understand what the issue is. Basically, the same covariance function written in two different forms, one will run and the other will not, even though they are exactly the same function. What am I missing? I paste the code for a simplified version of the model that reproduces this issue. Data are attached (n = 2). I am using cmdstanpy (stan version 2.23).

  functions{
        matrix matern52_v1(matrix qx,matrix qy, real s, real l){
            matrix[rows(qx),cols(qx)] r = sqrt(qx + qy)/l;
            return square(s)*(1 + sqrt(5)*r + 5/3*(r .* r)) .* exp(-sqrt(5)*r);
        }

        matrix matern52_v2(matrix qx,matrix qy, real s, real l){
            matrix[rows(qx),cols(qx)] r = sqrt(qx/square(l) + qy/square(l));
            return square(s)*(1 + sqrt(5)*r + 5/3*(r .* r)) .* exp(-sqrt(5)*r);
        }
    }
    data {
        int<lower=0> n;
        vector[n] lon;
        vector[n] lat;
        vector[n] y;
    }
    transformed data{
        matrix[n,n] qx;
        matrix[n,n] qy;
        
        for(i in 1:n){
            for(j in 1:n){
                qx[i,j] = (lon[i] - lon[j]) * (lon[i] - lon[j])/100;
                qy[i,j] = (lat[i] - lat[j]) * (lat[i] - lat[j])/100;
            }
        }
    }
    parameters {
        real<lower=0> l;
        real<lower=0> s;
    }
    model {
        matrix[n,n] q;
        matrix[n,n] L;
        l ~ std_normal();
        s ~ std_normal();
        q = matern52_v1(qx,qy,s,l);
        //q = matern52_v2(qx,qy,s,l);
        L = cholesky_decompose(q);
        y ~ multi_normal_cholesky(rep_vector(0,n),L);  
    }

Note that the functions ‘matern52_v1’ and ‘matern52_v2’ are exactly the same function just written differently. However, when I use matern52_v1 the code will run fine whereas when I use matern52_v2 the code will throw the following error:

RuntimeError: Error during sampling.
chain 1 returned error code 70

The *stdout.txt file contains the following message:

Rejecting initial value:
 Gradient evaluated at the initial value is not finite.
 Stan can't start sampling from this initial value.

lat.txt (50 Bytes)
lon.txt (52 Bytes)
y.txt (51 Bytes)

1 Like
2

I don’t have a solution, but I just wanted to let you know that I was able to recreate the error with cmdstan v2.27.0 using your code/data. So at the very least, it’s not some bizarre problem only on your machine.

2 Likes
3

Put a print statement inside both functions to see what r gets compute as, then call them on a small data input in the transformed data block.

4

The two functions return exactly the same values (i.e., up to 16 decimal places).

5

Hm. Ok, next step would be to run it on the full model but with the print still there. Make sure to use the same seed for both versions (or use both in the model block in a single model).