Rstan error in vb while using nngp

techniques

#1

Hi,

I am very new to Stan (Rstan). I was recently trying to model nngp using Lu’s example. I am giving here the full code as I added a weight factor while computing the covariance matrix. Initially I faced a problem on the new covariance matrix being “non-positive definite” in cholesky_decompose() after factoring in the weights. After I read few posts, I decided to force it to become pd and symmetric (which it already was). Although the problem is gone when I am running stan using NUTS (it has become slower though), I am getting an error while using the vb as:

Error in sampler$call_sampler(c(args, dotlist)) : Expecting a string vector: [type=NULL; required=STRSXP]

Can you please take a look and let me know how can I use variational bayes with advi here. I also don’t know if it is possible to make stan with nuts run faster. I am copying the R code as well. You can download the data from this link.

Thanks,
Subhomoy

 /* Latent NNGP model*/

  functions{
    
      matrix make_positive_definite(matrix m){
        
       int N = cols(m);
       matrix[N,N] es= eigenvectors_sym(m);
       vector[N] esv= eigenvalues_sym(m);
       real delta = 2 * 1e-5; 
       vector[N] tau; 
       matrix[N,N] dm;
       matrix[N,N] mdm;
       
         for(i in 1:N){
           tau[i] = fmax(0,delta-esv[i]);
         }
      dm = es * diag_matrix(tau) * es'; 
      
      return m + dm; 
      
     }

      
      real nngp_w1_lpdf(vector w1, real sigmasq1, real phi1, matrix NN_dist1,matrix NN_cor1,
                       matrix NN_distM1, matrix NN_corM1, int[,] NN_ind1, int N1, int M){

          vector[N1] V;
          vector[N1] I_Aw = w1;
          int dim;
          int h;

          for (i in 2:N1) {

              matrix[ i < (M + 1)? (i - 1) : M, i < (M + 1)? (i - 1): M]
              iNNdistM;
              matrix[ i < (M + 1)? (i - 1) : M, i < (M + 1)? (i - 1): M]
              iNNCholL;
              vector[ i < (M + 1)? (i - 1) : M] iNNcorr;
              vector[ i < (M + 1)? (i - 1) : M] v;
              row_vector[i < (M + 1)? (i - 1) : M] v2;

              dim = (i < (M + 1))? (i - 1) : M;

              if(dim == 1){iNNdistM[1, 1] = 1;}
              else{
                  h = 0;
                  for (j in 1:(dim - 1)){
                      for (k in (j + 1):dim){
                          h = h + 1;
                          iNNdistM[j, k] = exp(- phi1 * NN_distM1[(i - 1), h] * NN_corM1[(i - 1), h]);
                          iNNdistM[k, j] = iNNdistM[j, k];
                      }
                  }
                  for(j in 1:dim){
                      iNNdistM[j, j] = 1;
                  }
              }
              
              iNNdistM = 0.5*(iNNdistM + iNNdistM');
              iNNdistM = make_positive_definite(iNNdistM);

              iNNCholL = cholesky_decompose(iNNdistM);
              iNNcorr = to_vector(exp(- phi1 * (NN_dist1[(i - 1), 1:dim] .* NN_cor1[(i - 1), 1:dim]) ) );
              v = mdivide_left_tri_low(iNNCholL, iNNcorr);

              V[i] = 1 - dot_self(v);

              v2 = mdivide_right_tri_low(v', iNNCholL);

              I_Aw[i] = I_Aw[i] - v2 * w1[NN_ind1[(i - 1), 1:dim]];

          }
          V[1] = 1;
          return - 0.5 * ( 1 / sigmasq1 * dot_product(I_Aw, (I_Aw ./ V)) +
                          sum(log(V)) + N1 * log(sigmasq1));
      }
    
      
  }


  data {
      int<lower=1> N1;
      int<lower=1> N;
      int<lower=1> M;
      int<lower=1> P;
      vector[N] Y;
      matrix[N, P + 1] X;
      int NN_ind1[N1 - 1, M];
      matrix[N1 - 1, M] NN_dist1;
      matrix[N1 - 1, (M * (M - 1) / 2)] NN_distM1;
      matrix[N1 - 1, M] NN_cor1;
      matrix[N1 - 1, (M * (M - 1) / 2)] NN_corM1;
      matrix[N,N1] H1;
      vector[P + 1] uB;
      matrix[P + 1, P + 1] VB;
      real ss1;
      real st;
      real ap1;
      real bp1;
  }

  transformed data {
      cholesky_factor_cov[P + 1] L_VB;
      L_VB = cholesky_decompose(VB);
  }

  parameters{
      vector[P + 1] beta;
      real<lower = 0> sigma1;
      real<lower = 0> tau;
      real<lower = 0> phi1;
      vector<lower = 0>[N1] w1;
  }

  transformed parameters {
      real sigmasq1 = square(sigma1);
      real tausq = square(tau);
  }

  model{
      beta ~ multi_normal_cholesky(uB, L_VB);
      phi1 ~ gamma(ap1, bp1);
      sigma1 ~ normal(0, ss1);
      tau ~ normal(0, st);
      w1 ~ nngp_w1(sigmasq1, phi1, NN_dist1, NN_cor1, NN_distM1, NN_corM1, NN_ind1, N1, M);
      Y ~ normal(X * beta + H1 * w1, tau);
  }

# R 
rm(list = ls())
library(rstan)
#library(shinystan)
#library(spNNGP)       


#------------------------------ NNGP data load ---------------------------------#
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
load(file.choose())


#------------------------------ initialization vals & other vars ---------------#
myinits <-    list(list(beta = c(1, 2, 2), sigma1 = 10, tau = 10, phi1 = 100))
#               list(beta = c(2, 1, 3), sigma1 = 2, tau = 1, phi1 = .2),
#               list(beta = c(0.5, 1.5, .8), sigma1 = 3, tau = 0.5, phi1 = .5))

parameters <- c("beta", "sigmasq1", "tausq", "phi1")


#------------------------------ run stan Variational Bayes ---------------------#
stan_model <- stan_model(file = file.choose())

# ADVI
stan_fit_vb <- vb(stan_model, data=data_nngp, output_samples=1000, seed=123,init=myinits, pars=parameters)


#------------------------------ run MCMC with NUTS -----------------------------#

samples <- stan(
  file = "~/nngp_test.stan",
  data = data_nngp,
  init = myinits,
  pars = parameters,
  iter = 1000,
  chains = 1,
  thin = 100,
  cores=1,
  seed = 123
)