# Error transforming variable L: lub_free: Correlation variable is nan, but must be in the interval [-1, 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

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!

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.

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

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.