# L cholesky factor gets NA for fit

I have a STAN model that runs without divergences but that generates NA for the upper triangular matrix entries in the L // Cholesky factor of the correlation matrix (which are supposed to just be 0). This prevents bridgesampling from working for example.

any idea what it could be or what i am doing wrong?

``````data {

int<lower=1>  N;                   // Number of total observations
int<lower=1>  J[N];                // Subject-indicator per observation
int<lower=1>  K[N];                // Bandit-indicator per observation

// Data
int        Y[N];                   // Response (risky = 1, safe = 0)
real       R[N];                   // Outcome (win = 1, loss = 0)

}
transformed data {

int  NJ = max(J);                  // Number of total subjects
int  NK = max(K);                  // Number of total bandits

}
parameters {

// Subject parameters
vector[4]     theta_mu;            // Population-level effects
matrix[4,NJ]  theta_pr;            // Standardized subject-level effects

// Subject variances
cholesky_factor_corr[4] L;         // Cholesky factor of correlation matrix
vector<lower=0>[4] sigma;          // Subject-level standard deviations

}
transformed parameters {

vector[NJ]  b1;                    // Inverse temperature
vector[NJ]  a1;                    // Learning rate (positive)
vector[NJ]  a2;                    // Learning rate (negative)
vector[NJ]  q0;                    // Initial Q-value

// Construction block
{

// Rotate random effects
matrix[NJ,4] theta = transpose(diag_pre_multiply(sigma, L) * theta_pr);

// Construct random effects
b1 = (theta_mu[1] + theta[,1]) * 8;
a1 = Phi_approx(theta_mu[2] + theta[,2]);
a2 = Phi_approx(theta_mu[3] + theta[,3]);
q0 = Phi_approx(theta_mu[4] + theta[,4]);

}

}
model {
real delta;
real eta;

// Initialize Q-values
real  Q[NJ,NK] = to_array_2d(rep_matrix(q0, NK));

// Construct linear predictor
vector[N] mu;
for (n in 1:N) {

// Compute (weighted) difference in expected values
mu[n] = b1[J[n]] * (Q[J[n],K[n]] - 0.5);

// Compute prediction error
delta = R[n] - Q[J[n],K[n]];

// Assign trial-level learning rate
eta = R[n] * a1[J[n]] + (1-R[n]) * a2[J[n]];

// Update state-action values
Q[J[n],K[n]] += eta * delta;

}

// Likelihood
target += bernoulli_logit_lpmf(Y | mu);

// Priors
target += normal_lpdf(theta_mu | 0, 2);
target += std_normal_lpdf(to_vector(theta_pr));
target += student_t_lpdf(sigma | 3, 0, 1);
target += lkj_corr_cholesky_lpdf(L | 1);

}
generated quantities {

matrix[4,4] Omega = multiply_lower_tri_self_transpose(L);
vector[NJ]  ratio = (a1 - a2) ./ (a1 + a2);

}

``````