Covariance Matrix Not Symmetric for LKJ prior

I am attempting to deal with a short time series with four different data sets of surveys. Some data sets have repeated measures on the same respondents. I am using a normal regression model but when there are repeated measurements the likelihood is multivariate normal. I have reason to believe from data exploration that the variance is dependent on demographic characteristics. I want to put LKJ priors on the covariance matrices. When I run the stan code I get the following error:
“Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: multi_normal_lpdf: Covariance matrix is not symmetric. Covariance matrix[1,2] = 0, but Covariance matrix[2,1] = 2.9525 (in ‘modelc141bfc6bac_stantestpolicyv2’ at line 101)”
I suspect I have not specified the prior on the covariance correctly but I’m not sure what I’m doing wrong.

data {
  int<lower=0> N; // total sample size
  int<lower=1> D; // dimension of non-demographic covariates
  int<lower=1> K; // number of demographic covariates
  int<lower=1> T; //number of time points
  //group 1
  int<lower=1> N_g1; // sample size group 1
  vector[D] x_g1[N_g1]; // covariates group 1
  vector[K]  cvec_g1[N_g1]; // demographic covariates group 1
  real y_g1[N_g1]; // response group 1
  //group 2
  int<lower=1> N_g2; // sample size group 2
  matrix[D,2] x_g2[N_g2]; // covariates group 2
  vector[K]  cvec_g2[N_g2]; // demographic covariates group 2
  vector[2] y_g2[N_g2]; // response group 2
  //group 3
  int<lower=1> N_g3; // sample size group 3
  matrix[D,6] x_g3[N_g3]; // covariates group 3
  vector[K]  cvec_g3[N_g3]; // demographic covariates group 3
  vector[6] y_g3[N_g3]; // response group 3
  // group 4
  int<lower=1> N_g4; // sample size group 4
  vector[D] x_g4[N_g4]; // covariates group 4
  vector[K]  cvec_g4[N_g4]; // demographic covariates group 4
  real y_g4[N_g4]; // response group 4
  // vector of demographic group
  int<lower=1> C_g1[N_g1];
  int<lower=1> C_g2[N_g2];
  int<lower=1> C_g3[N_g3];
  int<lower=1> C_g4[N_g4];


// The parameters accepted by the model. Our model
parameters {
  vector[K] alpha; // incenterepts
  matrix [D,T] beta; // regresion coefficients
  vector<lower=0>[K] sigma_alpha; //variance of demographic coefficients
  vector<lower=0>[D] sigma_beta; //variance of non-demographic coefficients
  vector<lower=0>[32] sigma_error_s1; // variance of normal for s=1 for each covariance
  vector<lower=0>[32] sigma_error_s4;
  cholesky_factor_corr[2] group2cor[32]; // covariance of matrix for people observed 2
  cholesky_factor_corr[6] group3cor[32]; // covariance of matrix for people observed 6
  vector<lower=0>[2] sigmag2[32];
  vector<lower=0>[6] sigmag3[32];
  real<lower=0> group1var[32];
  real<lower=0> group4var[32];
  //real mu;
  //real<lower=0> sigma;
transformed parameters{
  //int my_subset[2] = {1,7};
  cholesky_factor_cov[2] covg2[32];
  cholesky_factor_cov[6] covg3[32];
  for(k in 1:32) covg2[k] = diag_pre_multiply(sigmag2[k],group2cor[k]);
  for(k in 1:32) covg3[k] = diag_pre_multiply(sigmag3[k],group3cor[k]);
// The model to be estimated. We model the output
// 'y' to be normally distributed with mean 'mu'
// and standard deviation 'sigma'.
model {
  //y ~ normal(mu, sigma);
  for(i in 1:K){
    sigma_alpha[i] ~ uniform(0,1);
  for(i in 1:D){
    sigma_beta[i] ~ uniform(0,1);
  for(i in 1:K){
    alpha[i] ~ normal(0, sigma_alpha[i]);
  for(i in 1:D){
    beta[i,:] ~ normal(0, sigma_beta[i]);
  for(i in 1:32){
    group2cor[i] ~ lkj_corr_cholesky(1);
    group3cor[i] ~ lkj_corr_cholesky(1);
    sigmag2[i] ~ uniform(0,1);
    sigmag3[i] ~ uniform(0,1);
    group1var[i] ~ uniform(0,1);
    group4var[i] ~ uniform(0,1);
  for(i in 1:N_g1){
    int c = C_g1[i];
    y_g1[i] ~ normal(dot_product(cvec_g1[i], alpha)+dot_product(x_g1[i],beta[:,1]), group1var[c]);
  for(i in 1:N_g2){
    int c = C_g2[i];
    vector[2] mu2;
    mu2[1] = dot_product(cvec_g2[i], alpha)+x_g2[i][:,1]'*beta[:, 1];
    mu2[2] = dot_product(cvec_g2[i], alpha)+x_g2[i][:,2]'*beta[:, 7];
    y_g2[i] ~ multi_normal(mu2, covg2[c]);
  for(i in 1:N_g3){
    int c = C_g3[i];
    vector[6] mu3;
    mu3[1] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,2]'*beta[:, 2];
    mu3[2] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,3]'*beta[:, 3];
    mu3[3] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,4]'*beta[:, 4];
    mu3[4] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,5]'*beta[:, 5];
    mu3[5] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,6]'*beta[:,6];
    mu3[6] = dot_product(cvec_g3[i], alpha)+x_g3[i][:,7]'*beta[:,7];
    y_g3[i] ~ multi_normal(mu3, covg3[c]);
  for(i in 1:N_g4){
    int c = C_g4[i];
    y_g4[i] ~ normal(dot_product(cvec_g4[i], alpha)+dot_product(x_g4[i],beta[:,1]), group1var[c]);


Your “covariance matrices” are Cholesky factors, not symmetric matrices.
Just use multi_normal_cholesky(mu, covg) instead of multi_normal(mu, covg).

If a parameter has uniform(0,1) prior then it should also have <lower=0,upper=1> bounds on its declaration. We usually don’t recommend uniform priors.