How to update one element of correlation each time

General
1

I have a correlation matrix with complex structure; for the estimation, I want to update one element in the correlation matrix each time because each element will have different values and each element’s updated values depends on the others. However, my current code assigns the same value to each updated element. My simplified code looks like this:

transformed parameters{
  corr_matrix[sum(K)] R;

  R = rep_matrix(0, sum(K), sum(K));
  for (m in 1:sum(K)){
    R[m, m] = 1;
  }
 
 // problems from here: this code assigns same value to all (m, n) elements, 
// but I want each element with different values from uniform (-1, 1);
  for (m in 1:K[1]){
    for (n in (K[1]+1): sum(K)){
        R[m, n] = nu; // want the nu to be different each time
        R[n, m] = R[m, n];
    }
  }
}

model {
  nu ~ uniform(-1, 1);
  print(R);
}
2

That code does not define a valid prior over the space of K \times K correlation matrices (if K > 2) because it puts positive probability on symmetric matrices that are not positive definite. Hence you will get a lot of rejections and NUTS will sample inefficiently.

3

The code I showed was a far much simplified version. My detailed method was
able to guarantee the correlation matrix to be positive definite. But I am just not sure how
to implement this in Stan properly with regards to update each element in the correlation matrix each time.

4

All the elements of the correlation matrix are changed at each leapfrog step. You just have to do your elementwise transformation correctly.

5

I wrote my code to estimate correlation matrix in R first and it worked well. Now I wanted to incorporate this part into the Stan code. I am not sure whether my replacement of runiform (0, 1) in the original R code with the “nu” parameter is proper. Also, the estimation of R from the following stan code is not as good as with my previous R code. Do you know what might caused the differences?

transformed parameters{
  corr_matrix[sum(K)] R;
  cov_matrix[sum(K)] VCOV;

  R = rep_matrix(0, sum(K), sum(K));
  for (m in 1:sum(K)){
    R[m, m] = 1;
  }

  // between-block correlations
  {
    matrix[sum(K), sum(K)] B_R;
    matrix[sum(K), sum(K)] R_N;
    matrix[sum(K), sum(K)] R_temp;
    matrix[sum(K), sum(K)] R_temp_C;
    matrix[sum(K), sum(K)] R_temp_N;

    real f_1;
    real f_minus1;
    real f_0;
    real a;
    real b;
    real c;
    real root1;
    real root2;
    real l_q1q2;
    real u_q1q2;
    real r_q1q2_N;
    real l_q1q21;
    real l_q1q22;
    real u_q1q21;
    real u_q1q22;
    real alpha_temp;
    real r_q1q2_C;
    real alpha_q1q2;
    real C_R;
    real ratio_N_C;
    real ratio_q1q2;

    R_N = R;
    B_R = rep_matrix(0, sum(K), sum(K));
    for (i in 1:N){
      B_R = B_R + inverse(diag_matrix(D)) * alpha[i] * (alpha[i]') * (inverse(diag_matrix(D))'); 
    }

    for (m in 1:K[1]){
      for (n in (K[1]+1): sum(K)){
        R_temp = R_N;
        R_temp[m, n] = 1;
        R_temp[n, m] = 1;
        f_1 = determinant(R_temp);
        R_temp[m, n] = -1;
        R_temp[n, m] = -1;
        f_minus1 = determinant(R_temp);
        R_temp[m, n] = 0;
        R_temp[n, m] = 0;
        f_0 = determinant(R_temp);
        a = (f_1 + f_minus1 - 2 * f_0) / 2; 
        b = (f_1 - f_minus1) / 2;
        c = f_0;

        if (a < 0){
          root1 = (-b + sqrt(b^2-4*a*c))/2/a;
          root2 = (-b - sqrt(b^2-4*a*c))/2/a;

          if (root1 > -1) l_q1q2 = root1; 
          else l_q1q2 = -1;

          l_q1q2 = ceil(l_q1q2*100)/100;

          if (root2 < 1) u_q1q2 = root2;
          else u_q1q2 = 1; 

          u_q1q2 = floor(u_q1q2*100)/100;
          r_q1q2_N = l_q1q2 + (u_q1q2 - l_q1q2) * nu; // nu is a random value between 0 and 1
        }


        if (a > 0){
          if ((b^2-4*a*c) < 0){
            r_q1q2_N = -1 + 2 * nu;  // nu is a random value between 0 and 1
          }
          if ((b^2-4*a*c) >= 0){
            root1 = (-b-sqrt(b^2-4*a*c))/2/a;
            root2 = (-b+sqrt(b^2-4*a*c))/2/a;
            l_q1q21 = -1;
            
            if (root1 > -1) u_q1q21 = root1; 
            else u_q1q21 = -1;

            u_q1q21 = floor(u_q1q21*100)/100;

            if (root2 < 1) l_q1q22 = root2; 
            else l_q1q22 = 1;

            l_q1q22 = ceil(l_q1q22*100)/100;
            u_q1q22 = 1;
            alpha_temp = nu; // nu is a random value between 0 and 1

            if (alpha_temp <= (u_q1q21-l_q1q21) / (u_q1q21-l_q1q21+u_q1q22-l_q1q22)) 
              r_q1q2_N = l_q1q21 + (u_q1q21 - l_q1q21 + u_q1q22-l_q1q22) * alpha_temp;
            if (alpha_temp > (u_q1q21-l_q1q21) / (u_q1q21-l_q1q21+u_q1q22-l_q1q22))
              r_q1q2_N=l_q1q22+(u_q1q21-l_q1q21+u_q1q22-l_q1q22)*alpha_temp-(u_q1q21-l_q1q21);

          }
        }

        R_temp_C = R_N;
        R_temp_N = R_N;
        R_temp_N[m, n] = r_q1q2_N;
        R_temp_N[n, m]= R_temp_N[m, n];
        r_q1q2_C = R_temp_C[m, n]; 
        alpha_q1q2 = nu; // nu is a random value between 0 and 1

        C_R = 1000^.5; 
        ratio_N_C = exp((log(determinant(R_temp_N)) - log(determinant(R_temp_C))) * (-N/2) - 1/2 * sum(diagonal(inverse(R_temp_N) * B_R)) + 1/2 * sum(diagonal(inverse(R_temp_C) * B_R))-(r_q1q2_N)^2/(C_R^2)/2+(r_q1q2_C)^2/(C_R^2)/2);
        
        if (ratio_N_C < 1) ratio_q1q2 = ratio_N_C; 
        else ratio_q1q2 = 1;

        if(ratio_q1q2 > alpha_q1q2){
          R_N[m, n] = r_q1q2_N;
          R_N[n, m] = R_N[m, n];
         } 

      }
    }

  R = R_N;
  }  

  VCOV = quad_form_diag(R, D);
}

model {
  nu ~ uniform(0, 1); // nu is the random value for the estimation of correlation matrix R
  for(n in 1:N){
    y[n] ~ multi_normal(zero_k, VCOV);
  }