Heteroskedastic noise ( again ! )

[Sébastien Coube]

Hello, 

I am coming back to you in order to ask (again !) how to model spatial heteroskedastic noise using INLA SPDE. The R code below shows what I mean; the objective is to estimate noise_variance and spatial_field. 

A few years ago, I have been in this forum to ask the same question. The answer that was kindly suggested to me was to run an INLA model with two Gaussian SPDEs : one for the Matérn field, one for the Gaussian noise field. The two SPDEs were differring in the size of their A-matrices, the noise SPDE having a diagonal A-matrix. The solution did work, but encountered difficulties when the number of observations and the number of variables explaining the noise variance augmented. 

I would like to try again, and to be sure that I am not missing a new method that has been introduced in INLA since the last time. 

Again, thank you all for the terrific work you are doing ! 

Kind regards

Sébastien Coube-Sisqueille

# spatial locations
locs = as.matrix(expand.grid(seq(1e2), seq(1e2)))

# Gaussian field explaining the mean
spatial_field = matrix(rnorm(1e4), 1e2)
sqexp_cov_chol = exp(-as.matrix((dist(seq(1e2), diag = T)/10)^2))
diag(sqexp_cov_chol) = 1.001
sqexp_cov_chol = chol(sqexp_cov_chol)
spatial_field = t(sqexp_cov_chol) %*% spatial_field %*% sqexp_cov_chol
plot(
  locs,
  col = heat.colors(100)[cut(spatial_field, seq(min(spatial_field), max(spatial_field), length.out = 101))],
  pch = 15,
  main = "spatial field"
     )

# Gaussian field explaining the mean
noise_variance = matrix(rnorm(1e4), 1e2)
noise_variance = - 2 + t(sqexp_cov_chol) %*% noise_variance %*% sqexp_cov_chol
noise_variance = exp(noise_variance)
noise_realization = sqrt(noise_variance) * rnorm(1e4)
plot(
  locs,
  col = heat.colors(100)[cut(noise_variance, seq(min(noise_variance), max(noise_variance), length.out = 101))],
  pch = 15,
  main = "variance of the noise"
     )
plot(
  locs,
  col = heat.colors(100)[cut(noise_realization, seq(min(noise_realization), max(noise_realization), length.out = 101))],
  pch = 15,
  main = "realization of the noise"
     )

# combining spatial field and noise
observed_field = spatial_field + noise_realization
plot(
  locs,
  col = heat.colors(100)[cut(observed_field, seq(min(observed_field), max(observed_field), length.out = 101))],
  pch = 15,
  main = "observed field"
     )