Hello,

I am trying to generate a variance-covariance matrix using the decomposition from Barnard et al 2000, where the variance-covariance matrix can be written as a product of a diagonal matrix and a correlation matrix. I want to draw the off-diagonal values of the correlation matrix from a (shifted and scaled) Beta distribution where the shape parameters are not necessarily equal to each other.

The problem I am running into is ensuring that the correlation matrix is positive definite. Iām wondering if there is a way to implement this restriction in Stan. I keep getting this error during the sampling process:

"Chain 1: Rejecting initial value:

Chain 1: Error evaluating the log probability at the initial value.

Exception: validate transformed params: S[i_0__] is not positive definite.ā

Here is what my code currently looks like (I omitted some of the lines that arenāt directly related to this problem):

```
// construct a correlation matrix
functions {
matrix to_corr(vector x, int K) {
matrix[K, K] y;
int pos = 1;
for (i in 1: K) {
for (j in 1:(i - 1)) {
y[i, j] = x[pos];
y[j, i] = x[pos];
pos += 1;
}
y[i, i] = 1;
}
//print(y);
return y;
}
}
data {
int<lower=1> Q; // dimension of the v-cov matrices
int<lower=1> I; // no. of subjects
}
transformed data {
int<lower=0> off_diag_Q = (Q*(Q-1))/2;
}
parameters {
vector[Q] log_Q_Si[I];
vector <lower=0, upper=1>[off_diag_Q] r_raw[I]; //values to draw from Beta distribution
real<lower=0> hyper_alpha; // hyperparameter on the Beta prior
real<lower=0> hyper_beta; //hyperparameter on the Beta prior
}
transformed parameters{
cov_matrix[Q] S[I]; // subject specific vcovās
vector<lower=0>[Q] D_Si[I]; //non-negative diagonals of the Si decomposition
matrix[Q,Q] R[I]; // corr matrix of Si decomposition
for(i in 1:I){
R[i] = to_corr(2 * r_raw[i] - 1, Q)} //shift & scale to get valid correlations between -1 and 1
D_Si[i] = exp(log_Q_Si[i]);
S[i] = quad_form_diag(R[i], D_Si[i]); //get back the S_i's
}
}
model {
//hyperparameters on the S_is
hyper_alpha ~ exponential(0.1);
hyper_beta ~ exponential(0.1);
for(i in 1:I) {
r_raw[i] ~ beta(hyper_alpha, hyper_beta);
}
}
```

(If it helps, I only need this for the 3x3 matrix case.)

Thank you!