Simple spatial guassian process model - unable to recover parameters

Modeling
1

Hello,

I am tryin to fit a simple spatial GP model - first simulating the field with the R package geoR , function geoR::grf using an exponential covariance function with variance and range parameters (0.8 and 0.5) respectively. Then fitting the model in Stan does not recover the parameters (nor their product) - so I am wondering if I am doing something wrong.

Data simulated in R:

library(geoR)
library(rstan)
pref_field <- geoR::grf(
  200,
  grid = "reg",
  nsim = 1,
  cov.model = "exponential",
  cov.pars = c(0.8,0.5),
  kappa=0,
  message=FALSE
) 

stan_data <- list(
  s= length(pref_field$data),
  coords = pref_field$coords,
  phi_locs = pref_field$data
)

Stan model:

data {
  int<lower=0> s;
  array[s] vector[2] coords;
  vector[s] phi_locs;
}
parameters {
  real<lower=0> sigma2;
  real<lower=0> rho;
  real phi_mean;
}
transformed parameters{
  matrix[s,s] SIGMA;
  matrix[s,s] L_SIGMA;
  
  SIGMA = gp_exponential_cov(coords, sqrt(sigma2), rho);
  L_SIGMA = cholesky_decompose(SIGMA);
  
}
model {
  target += normal_lpdf(phi_mean|0,100);
  target += gamma_lpdf(sigma2 | 2, 0.25);
  target += inv_gamma_lpdf(rho | 3, 1);
  target += multi_normal_cholesky_lpdf(phi_locs|rep_vector(phi_mean,s),L_SIGMA);
}

Results:

> round(summary(simple_pp,c("sigma2","rho","phi_mean"))$summary[,c(1,4,7)],3)
                mean          2.5%      75%
sigma2       1.758          0.466  2.030
rho              1.051         0.263  1.199
phi_mean   -0.246        -2.563  0.293

Thanks for any advice!

2

The Stan code generally looks correct, though using the phi_mean ~ normal(0,100) sampling notation is probably preferred. It looks like the posterior credible intervals do contain the true parameter values. I would first suggest retrying with more simulated data. It’s not necessarily guaranteed that the posterior means will exactly recover the generative parameters (only in expectation). If problems persist here are some suggestions for troubleshooting:

  • Normal(0,100) is a very wide prior for the spatial process mean, implying that it could plausibly be between about -200 and 200. Does this make sense for your simulated data? Does the problem go away if you use e.g. Normal(0,1)?
  • Try directly simulating the field without relying on geoR or double checking with a different package to verify there isn’t an issue with that package or a mistranslation of the process parameters
2 Likes
3

Depends who you ask! I like the distribution notation. To help, we renamed it from “sampling statement” to “distribution statement” because it’s normally read “is distributed as”. Many folks, including some of our developers, are lobbying us to remove the ~ notation altogether because they find it’s confusing users because it’s not actually generating phi_mean by sampling a normal(0, 100). The exact equivalent to the sampling statement with target += is

target += normal_lupdf(phi_mean | 0, 100);

where the lupdf has the extra u for unnormalized, which drops constants like the distribution statement.

Like @js592, I’d also recommend tighter priors. The gamma prior on sigma2 has a mean of 8 (2 / 0.25) and a standard deviation of roughly 6 in case you thought we were using a different parameterization of gamma than shape and rate (aka inverse scale).

Otherwise, I’d recommend @js592’s advice of doing the simulation directly or in Stan. At the lowest level, I’d at least check that the different distributions have the same parameterizations.

1 Like
4

Ah - indeed I mistakenly thought that the gamma distribution was parametrized with shape and scale. Thanks for catching that.
Thanks both for the recommendations. That helped.