Location parameter Inf in Latent Gaussian and long warmup

Modeling
1

Hi,

I am trying to fit a latent Gaussian model where the main objective is to obtain the posterior latent vector sf (see below in my code). In the toy set up, it is a very simple model where covariance matrices are diagonal; data and latent vectors are all normal. However, although the following stan code compiles okay, I am getting the warning "Exception: multi_normal_cholesky_lpdf: Location parameter[1] is inf, but must be finite! " while running the sampler. I think because of this the warmup stage is taking a huge time to go through. I am not that familiar with stan. It might be that the model set up is proper. I was hoping if anyone can help!

Below here is the code, and R code. I am running it in R with cmdstanr. Here is the dataset link.

Thanks,
SG

 /*Latent Gaussian*/

  functions{
     matrix cov_custom_LQ(real tausq, real sigmasq, vector sill, int p){

          vector[p] sillk= rep_vector(tausq,p) + sigmasq*sill;
          matrix[p, p] Q = diag_matrix(sillk);
          return cholesky_decompose(Q);          
      }
      
    matrix cov_custom_LR(real sigmasqR, vector sillR, int n){

          vector[n] sillk= sigmasqR*sillR;
          matrix[n, n] R = diag_matrix(sillk);
          return cholesky_decompose(R);          
      }
    }

  data {
      int<lower=1> n;
      int<lower=1> p;
      vector[n] y;
      matrix[n,1] HX;
      matrix[p,1] X;
      matrix[n,p] H;
      vector[p] sillQ;
      vector[n] sillR;
  }

   transformed data {
     vector[p] mu= rep_vector(0,p);
  }

  parameters{
      vector<lower = 0 >[1] beta;
      real<lower=0> sigmasq;
      real<lower=0> sigmasqR;
      real<lower=0> tausq;
      vector[p] s;
  }
  
  

  model{
      matrix[p,p] QL = cov_custom_LQ(tausq,sigmasq,sillQ,p);
      matrix[n,n] RL = cov_custom_LR(sigmasqR,sillR,n);
      beta ~ normal(1,5);
      sigmasqR ~ cauchy(0,1);
      tausq ~ cauchy(0,1);
      sigmasq ~ cauchy(0,1);
      s ~ multi_normal_cholesky(mu,QL);
      y ~ multi_normal_cholesky(HX*beta + H*s,RL);
  }
  
generated quantities {
   vector[p] sf;
   sf = s + X*beta;
   vector[n] ypost;
   ypost = H*sf;
}

Running code:

rm(list = ls())
library(cmdstanr)


#------------------------------ data load ---------------------------------#
load(file.choose())

#------------------------------ initialization vals & other vars ---------------#
myinits <-    list(list(beta = c(1), sigmasq = 1, tausq = 1, sigmasqR=1,s=rep(1,156)))
parameters <- c("beta", "sigmasq", "tausq","sigmasqR","s")

#------------------------------ run  -----------------------------#

file <- file.path("")
mod <- cmdstan_model(file)


fit <- mod$sample(
  data = data,
  seed = 123,
  chains = 1,
  parallel_chains = 1,
  refresh = 500,
  init = myinits
)
2

Is the model producing warnings about divergences as well?

3

I actually haven’t gone past the warm-up step as it was taking too much time. It produced the same warning many times. However, I think, I am able to bypass this problem by using non-centered parameterization on the latent vector’s distribution. I also started with a different initial value of eta (see below). However, I am not sure why the problem occurred. It would be good to figure out the reasons because, in the actual problem, my covariance matrices are more complex. Below here is the new version:

 /*Latent Gaussian*/

  functions{
     matrix cov_custom_LQ(real tausq, real sigmasq, vector sill, int p){

          vector[p] sillk= rep_vector(tausq,p) + sigmasq*sill;
          matrix[p, p] Q = diag_matrix(sillk);
          return cholesky_decompose(Q);          
      }
      
    matrix cov_custom_LR(real sigmasqR, vector sillR, int n){

          vector[n] sillk= sigmasqR*sillR;
          matrix[n, n] R = diag_matrix(sillk);
          return cholesky_decompose(R);          
      }
    }

  data {
      int<lower=1> n;
      int<lower=1> p;
      vector[n] y;
      matrix[n,1] HX;
      matrix[p,1] X;
      matrix[n,p] H;
      vector[p] sillQ;
      vector[n] sillR;
  }

   transformed data {
     vector[p] mu= rep_vector(0,p);
  }

  parameters{
      vector<lower = 0 >[1] beta;
      real<lower=0> sigmasq;
      real<lower=0> sigmasqR;
      real<lower=0> tausq;
      vector[p] eta;
  }
  
  transformed parameters {
    matrix[p,p] QL = cov_custom_LQ(tausq,sigmasq,sillQ,p);
    matrix[n,n] RL = cov_custom_LR(sigmasqR,sillR,n);
    vector[p] s= X*beta + QL*eta;
    vector[n] ymu=H*s;
  }
  

  model{
      beta ~ normal(1,5);
      sigmasqR ~ normal(0,10);
      tausq ~ normal(0,10);
      sigmasq ~ normal(0,1);
      //s ~ multi_normal_cholesky(mu,QL);
      eta ~ std_normal();
      y ~ multi_normal_cholesky(ymu,RL);
  }