Hi Stan Community,
Sorry for the long post – but at least the model is short!
I was hoping for some help learning how to use Stan to develop a hierarchical spatial model that infers a complete spatial field from noisy, biased, and gappy observations. Specifically, I’m trying to recreate the simple example here (Section 4). There is even a Matlab implementation of a custom Gibbs Sampler (zip file here;/GHCN_CONUS_Demo; note the link in the previous document does not work).
I’ve tried to start very basic, by regarding the posterior medians of the parameters and the process (from the Matlab solution) as ‘truth’, giving them as data to Stan, and then just letting it infer the process. This produces the correct result and is reasonably fast.
However, as soon as I start letting it infer the parameters it takes quite awhile, despite all the priors being conjugate. For example, the attached model took 1.5 hours for 500 iterations. The Gibbs sampler in Matlab takes about a minute for 1000 iterations.
Surely there must be something wrong with how I’m specifying the model. I suppose there could also be a problem with this debugging strategy, but I also have tried to specify the model completely with the real input data, as well as noise added to the ‘true’ posterior process, and have the same problems. I’ve also tried specifying an upper bound on the covariance parameters but this didn’t make a difference.
The model code I’m attaching is a demonstration of inferring just the sig_2 parameter. As noted, it took 1.5 hours for 500 iterations, produced a reasonable looking posterior process, but substantially underestimated the parameter : the median and max of sig_2 are 7.9 with a max of 9.8, while the posterior median ‘truth’ is 16.7.
A specific question, that perhaps betrays my lack of understanding of Stan, is: why am I getting this warning:
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is 0.
I don’t see how sig_2 * R could not be positive definite. This only happens early in the simulation, I guess it could be drawing a 0, but still seems strange.
Another problem I have is that the model hangs after finishing sampling. This happens when I specify the full covariance matrix (SIG_2 in code) in the transformed parameter. My guess is that it is struggling with the size of the matrix (730 x 730 x 500), but I don’t think this is THAT big…
Any general help would be much appreciated, including pointing to any examples/papers implementing such a model in Stan, which I haven’t been able to find. I can provide the pystan driver and data if you’d like, as well as figures comparing the posterior between Stan and the Gibbs Sampler in Matlab.
data {
int<lower=0> N; // process locations to infer at
vector[N] z; // temperature observations
real mu; // process mean
real<lower=0.0> tau_2; // data covariance
cov_matrix[N] R; // correlation matrix (scaled distances between locations)
}
transformed data {
vector[N] MU;
// cov_matrix[N] TAU_2;
MU = rep_vector(mu, N);
// TAU_2 = diag_matrix(rep_vector(tau_2, N));
}
parameters {
vector[N] y; // process
real<lower=0> sig_2; // process covariance
}
transformed parameters {
//cov_matrix[N] SIG_2;
//SIG_2 = sig_2 * R;
}
model {
// priors
// mu ~ normal(0, 25);
// tau_2 ~ inv_gamma(5, 25);
sig_2 ~ inv_gamma(5, 100);
// process/state level
y ~ multi_normal(MU, sig_2*R);
// data/observation level
//z ~ multi_normal(y, TAU_2); // believe this produces same result as below
z ~ normal(y, sqrt(tau_2);
}
Thanks,
CC