Testing…

Can you post some visuals of your raw data? Just y as a function of x with geom_point

That could help, I see. But it’s a bit difficult to share publicly the actual data or plot of the actual data of an ongoing research. Makes sense?

Ok, then maybe step through chunks of f with the pairs plot (you can ask for subsets of f via `c(f[1], f[2], …)`

). Maybe also add some posterior predictive checking, which might highlight particular data points for which the model is misspecified.

Thanks. I will do that. Thanks for the input. It was very helpful.

@andrjohns, any thoughts how to solve this?

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for:

**normal_rng(matrix, real)**

Available argument signatures for normal_rng:

normal_rng(real, real)

normal_rng(real, real[ ])

normal_rng(real, vector)

normal_rng(real, row_vector)

normal_rng(real, int)

normal_rng(real, int[ ])

normal_rng(real[ ], real)

normal_rng(real[ ], real[ ])

normal_rng(real[ ], vector)

normal_rng(real[ ], row_vector)

normal_rng(real[ ], int)

normal_rng(real[ ], int[ ])

normal_rng(vector, real)

normal_rng(vector, real[ ])

normal_rng(vector, vector)

normal_rng(vector, row_vector)

normal_rng(vector, int)

normal_rng(vector, int[ ])

normal_rng(row_vector, real)

normal_rng(row_vector, real[ ])

normal_rng(row_vector, vector)

normal_rng(row_vector, row_vector)

normal_rng(row_vector, int)

normal_rng(row_vector, int[ ])

normal_rng(int, real)

normal_rng(int, real[ ])

normal_rng(int, vector)

normal_rng(int, row_vector)

normal_rng(int, int)

normal_rng(int, int[ ])

normal_rng(int[ ], real)

normal_rng(int[ ], real[ ])

normal_rng(int[ ], vector)

normal_rng(int[ ], row_vector)

normal_rng(int[ ], int)

normal_rng(int[ ], int[ ])

error in ‘modelae463fb7c6fnew’ at line 47, column 52

```
47: y_rep[i] = normal_rng(a + L_K*eta[i], sigma);
^
```

The quick fix is to put parentheses around the computation of the latent GP:

```
y_rep[i] = normal_rng((a + L_K*eta)[i], sigma);`
```

But that’s fairly compute-inefficient, as it means computing the `(a + L_K*eta)`

repeatedly. Best to compute

```
vector[N] f = a + L_K*eta ;
vector[N] y_rep ;
for( i in 1:N ){
y_rep[i] = normal_rng( f[i] , sigma ) ;
}
```

Try using the vectorised rng instad:

```
generated quantities {
real y_rep[N] = normal_rng(a + L_K*eta, sigma);
}
```

Even better, you can integrate out `eta`

and `a`

analytically and go from N+4 dimensional posterior to 3 dimensional posterior. Specially, integrating out `eta`

is likely to remove the divergences. Model code looks like

```
model {
matrix[N, N] K_f = gp_exp_quad_cov(xn, sigma_f, lengthscale_f) +
sigma_intercept;
matrix[N, N] L_f = cholesky_decompose(add_diag(K_f, sigma));
yn ~ multi_normal_cholesky(zeros, L_f);
lengthscale_f ~ normal(0, 1);
sigma_f ~ normal(0, 1);
sigman ~ normal(0, 1);
}
```

where `sigma_intercept`

and `zeros`

are defined in transformed data

```
real sigma_intercept = 1;
vector[N] zeros = rep_vector(0, N);
```

The full code can be found in a case study Gaussian process demonstration with Stan

Very cool! Can’t believe I haven’t seen this trick until now; I clearly wasn’t paying sufficient attention last time I read your case study.

Indeed. Divergences are gone now. Thanks.

This is how I learned Gaussian processes. Put everything in the covariance matrix that you can. Integrated Laplace approximation will make it possible to use this approach also for non-normal observation models.

Excellent! Great that your covariance matrices are small enough for sufficiently fast computation.