Multi-block correlation matrix

Modeling
7

Thank you, Caesoma! I figured out a way to estimate this special correlation structure, and I was able to implement it in R. 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 problem?

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