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