Help reparameterize GP model to remove divergent transitions

Modeling
, , ,
21

Testing…

22

sigma ~ weibull(2,1) did not remove divergent transitions. Here’s the new mcmc_pairs. Any thoughts?

image
image639×860 31 KB

23

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

24

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?

25

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.

1 Like
26

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

27

@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);
                         ^
28

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 ) ;
}
1 Like
29

Try using the vectorised rng instad:

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

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

3 Likes
31

Thanks @avehtari. This is a great help. Let me apply it and hopefully remove the divergences.

32

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.

1 Like
33

Indeed. Divergences are gone now. Thanks.

3 Likes
34

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.

4 Likes