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