# L cholesky factor gets NA for fit

**URL:** https://discourse.mc-stan.org/t/l-cholesky-factor-gets-na-for-fit/29946
**Category:** Modeling
**Tags:** fitting-issues, rstan
**Created:** [January 5, 2023, 9:19am UTC](https://discourse.mc-stan.org/t/l-cholesky-factor-gets-na-for-fit/29946 "2023-01-05T09:19:24Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![levis](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/levis/32/8315_2.png) [@levis](https://discourse.mc-stan.org/u/levis)
#### Post date: [January 5, 2023, 9:19am UTC](https://discourse.mc-stan.org/t/l-cholesky-factor-gets-na-for-fit/29946/1 "2023-01-05T09:19:24Z")

</div>

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?

```plaintext
data {

    // Metadata
    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);
    
}

```