Error transforming variable L: lub_free: Correlation variable is nan, but must be in the interval [-1, 1]

General
1

I am trying to sample from the posterior of Gaussian copula model with Bernoulli marginals, the dimension of the multivariate Gaussian in this model is 542, and I want to sample from the posterior of the correlation matrices.

I have the following error messages:

    SAMPLING FOR MODEL '92fbbf86083410150a497cd7bdd2db93' NOW (CHAIN 1).
    Chain 1: Unrecoverable error evaluating the log probability at the initial value.
    Chain 1: Error transforming variable L: lub_free: Correlation variable is nan, but must be in the 
    interval [-1, 1]

The stan model is:

   functions{
     vector [] fFindBounds(row_vector vU, matrix L, int[] vY, vector vProb){
   int K = cols(vU);
   vector[K] mBounds;
   vector[K] lLogProbU;
   vector[K] u_star_lower;
   vector[K] u_star_upper;
   // real z_star_lower;
   real z_star;
   real z_star_L;
   real z_star_L_sign;
   vector[K] z;
   vector[K] v;
   vector[K] out[2];
   for(k in 1:K){
   int km1 = k - 1;

   // upper bound only
   mBounds[k] = inv_Phi(vProb[k]);
   z_star_L = (mBounds[k] - ((k > 1) ? L[k,1:km1] * head(z,km1) : 0));
   z_star_L_sign = (z_star_L < 0) ? -1 : 1 ;
   if (is_inf(z_star_L) || is_nan(z_star_L)){
   z_star_L = 1e8;
   };
   // print(z_star_L_sign, \" \", z_star_L, \" \",L[k,k]);
           z_star = (z_star_L * z_star_L_sign) / ((L[k,k] < 1e-8) ? 1e-8 : L[k,k]);
         if(vY[k] == 0){
           u_star_lower[k] = 0;
           u_star_upper[k] = Phi(z_star);
         }else{ 
           // lower bound only
           u_star_lower[k] = Phi(z_star);
           u_star_upper[k] = 1;
         }
         v[k] = u_star_lower[k] + (u_star_upper[k]-u_star_lower[k]) * vU[k];
         z[k] = inv_Phi(v[k]);
         lLogProbU[k] = u_star_upper[k] - u_star_lower[k];
       }

       for (k in 1:K){
       lLogProbU[k] = (lLogProbU[k] < 1e-6) ? 1e-6 : lLogProbU[k];
       }
       
       out[1] = lLogProbU;
   //    print(log(lLogProbU));
       out[2] = v;
       
       return(out);
     }
   }


   data{
   int N;
   int K;
   int Y[N,K];
   }

   //transformed data{
   //vector[K] vZeros;
   //for(i in 1:K)
   //vZeros[i] = 0;
   //}

   parameters{
   cholesky_factor_corr[K] L; // correlation between Xs
   vector<lower=1e-6,upper=1 - 1e-6>[K] probs; // probability of 0 of marginal bernolli distributions


   // augmented likelihood parameters
   matrix<lower=1e-6,upper=1-1e-6>[N,K] u;
   }


   transformed parameters{
   //matrix[N,K] x;
   //vector[N] lLogProb;
   matrix[N,K] mLogProbU;
   //matrix[N,K] z;
   //vector[K] lOut[N,2];
   //matrix[N,K] mv;




   for(i in 1:N){
   mLogProbU[i] = to_row_vector(fFindBounds(u[i],L,Y[i],probs)[1]);
   }

   //for(i in 1:N){
   //lLogProb[i] = multi_normal_cholesky_lpdf(to_vector(x[i])|vZeros,L);
   //}
   }

   model{
   L ~ lkj_corr_cholesky(1.0);

   for(i in 1:N){
   for(j in 1:K)
   target += log(mLogProbU[i,j]);
   print(target());
   }


   // probs ~ normal(5,10);
   }

   generated quantities{
   corr_matrix[K] Sigma;
   matrix[N,K] x;
   matrix[N,K] z;
   matrix[N,K] mv;



   for(i in 1:N){
   mv[i] = to_row_vector(fFindBounds(u[i],L,Y[i],probs)[2]);
       for(j in 1:K){
   z[i,j] = inv_Phi(mv[i,j]);
   }
   }

   for(i in 1:N){
   x[i] =  to_row_vector(L * to_vector(z[i]));
   }

   Sigma = L * L';
   }

The r codes for sampling:

    biPCAmodel <- stan_model(model_code = model_discrete_Gaussian_copula_family)
    # data preprocessing
    library(readxl)
    senators<-read_xls("senate_voting_data.xls")

    senators.names<-names(senators)[-c(1,2)]
    rev.party.state.names<-lapply(X=strsplit(gsub(pattern="[.]",replacement="",x=senators.names),split=" "),FUN = rev)

    senators.party <- lapply(X = rev.party.state.names, FUN = function(x)(unlist(x)[1]))
    senators.party <- unlist(senators.party)

    senators.last.names <- lapply(X = rev.party.state.names, FUN = function(x)(unlist(x)[4]))
    senators.last.names <- unlist(senators.last.names)

    senators.last.names.party <- paste(senators.last.names, senators.party, sep = '-')
    senators_new <- as.data.frame(t(senators[,-c(1,2)]))

    colnames(senators_new) <- NULL
    rownames(senators_new) <- NULL

    senators_bin <- senators_new
    senators_bin[as.matrix(senators_new) == -1] <- 0



    probs_init <- apply(senators_bin, MARGIN = 2, function(x) sum(x)/length(x))

    res <- sampling(biPCAmodel, list(N = nrow(senators_bin), K = ncol(senators_bin), Y = senators_bin), verbose = T, pars = c("Sigma", "x", "probs"), init = list(
      list(probs = probs_init), list(probs = probs_init), list(probs = probs_init), list(probs = probs_init) 
    ))

It is like this is the problem of the high dimension of correlation matrix. I’ve also added some checking in the fFindBounds function to avoid producing inf and nan, but I cannot figure out why there is error with the Cholesky factor martix L. Can anyone help to find any mistakes in the model? Thanks!

2

When something is declared as a Cholesky factor in the parameters block and you get messages like this, it is usually the case that the gradient evaluated to NaN due to numerical problems.

3

Thanks Ben. So is there any way to avoid this? I cannot get any samples now because of this error.

4

I would first try setting init_r to some value less than its default of 2 but ultimately you may have to write a more numerically stable function.