Ok, so I flipped the problem a bit and I think this works. Let’s get a PD matrix that is as close to possible to the correlations you want.

```
functions{
matrix angle_vec2mat (vector angle, int K) {
int N = num_elements(angle);
matrix[K, K] mat = add_diag(identity_matrix(K), rep_vector(-1, K));
int count = K;
int pos = 1;
for (k in 1:K - 1){
count -= 1;
mat[k + 1:K, k] = segment(angle, pos, count);
pos += count;
}
return mat;
}
matrix angle2cor (vector angle, int N){
matrix[N, N] angle_mat = angle_vec2mat(angle, N);
matrix[N, N] inv_mat = add_diag(angle_mat, rep_vector(1, N));
inv_mat[, 1] = cos(angle_mat[, 1]);
for (i in 2:N) {
inv_mat[i, i] = prod(sin(angle_mat[i, 1:(i - 1)]));
if (i > 2) {
for (j in 2:i - 1 ) {
inv_mat[i, j] = cos(angle_mat[i, j]) * prod(sin(angle_mat[i, 1:(j - 1)]));
}
}
}
return inv_mat;
}
vector lower_elements(matrix M, int tri_size){
int n = rows(M);
int counter = 1;
vector[tri_size] out;
for (i in 2:n){
for (j in 1:(i - 1)) {
out[counter] = M[i, j];
counter += 1;
}
}
return out;
}
}
data{
int d;
int n;
vector[d] dat[n];
real corpriorscale;
}
transformed data {
int M = (d * d - d) / 2;
}
parameters{
vector<lower=0, upper=pi()>[M] rawcor;
vector<lower=-1, upper=1>[M] cors;
vector[d] logsdvec;
vector[d] mu;
}
transformed parameters{
cholesky_factor_corr[d] L_m = angle2cor(rawcor,d);
corr_matrix[d] mcor = multiply_lower_tri_self_transpose(L_m);
}
model{
cors ~ normal(lower_elements(mcor, M), corpriorscale);
logsdvec ~ normal(0,10);
mu ~ normal(0, 10);
dat ~ multi_normal_cholesky(mu, diag_pre_multiply(exp(logsdvec), L_m));
}
```

We can really stress test this by flipping the `cors`

parameters into data. Just like you were doing when you were exposing functions. Though we can’t do that anymore because we need the soft constraint of the normal to regularize a PD near the ones we’re looking for.

```
functions{
matrix angle_vec2mat (vector angle, int K) {
int N = num_elements(angle);
matrix[K, K] mat = add_diag(identity_matrix(K), rep_vector(-1, K));
int count = K;
int pos = 1;
for (k in 1:K - 1){
count -= 1;
mat[k + 1:K, k] = segment(angle, pos, count);
pos += count;
}
return mat;
}
matrix angle2cor (vector angle, int N){
matrix[N, N] angle_mat = angle_vec2mat(angle, N);
matrix[N, N] inv_mat = add_diag(angle_mat, rep_vector(1, N));
inv_mat[, 1] = cos(angle_mat[, 1]);
for (i in 2:N) {
inv_mat[i, i] = prod(sin(angle_mat[i, 1:(i - 1)]));
if (i > 2) {
for (j in 2:i - 1 ) {
inv_mat[i, j] = cos(angle_mat[i, j]) * prod(sin(angle_mat[i, 1:(j - 1)]));
}
}
}
return inv_mat;
}
vector lower_elements(matrix M, int tri_size){
int n = rows(M);
int counter = 1;
vector[tri_size] out;
for (i in 2:n){
for (j in 1:(i - 1)) {
out[counter] = M[i, j];
counter += 1;
}
}
return out;
}
}
data{
int d;
int n;
vector[d] dat[n];
real corpriorscale;
real x;
}
transformed data {
int M = (d * d - d) / 2;
vector[M] cors = rep_vector(x, M);
}
parameters{
vector<lower=0, upper=pi()>[M] rawcor;
// vector<lower=-1, upper=1>[M] cors;
vector[d] logsdvec;
vector[d] mu;
}
transformed parameters{
cholesky_factor_corr[d] L_m = angle2cor(rawcor,d);
corr_matrix[d] mcor = multiply_lower_tri_self_transpose(L_m);
}
model{
cors ~ normal(lower_elements(mcor, M), corpriorscale);
logsdvec ~ normal(0,10);
mu ~ normal(0, 10);
dat ~ multi_normal_cholesky(mu, diag_pre_multiply(exp(logsdvec), L_m));
}
```

Let’s test the original where `corpriorscale = 10`

and the correlation is 0 so `x = 0`

. Then where `corpriorscale = 0.2`

and `x = 0.9`

to show how it pushes the correlation matrix up toward 0.9.

```
sdat <- list(n=n,d=d, dat=dat,corpriorscale= 10,
x = 0)
> mod_out$summary("mcor")
# A tibble: 25 x 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mcor[1,1] 1 1 0 0 1 1 NA NA NA
2 mcor[2,1] 0.722 0.734 0.0907 0.0874 0.551 0.844 1.00 2133. 2383.
3 mcor[3,1] -0.172 -0.179 0.177 0.178 -0.453 0.130 1.00 4418. 2912.
4 mcor[4,1] 0.239 0.249 0.173 0.172 -0.0664 0.507 1.00 2147. 2321.
5 mcor[5,1] -0.409 -0.422 0.153 0.151 -0.641 -0.131 1.00 3097. 2750.
6 mcor[1,2] 0.722 0.734 0.0907 0.0874 0.551 0.844 1.00 2133. 2383.
7 mcor[2,2] 1 1 0 0 1 1 NA NA NA
8 mcor[3,2] 0.0120 0.0155 0.182 0.190 -0.293 0.303 1.00 4187. 2889.
9 mcor[4,2] 0.828 0.837 0.0600 0.0536 0.718 0.907 1.00 3130. 2966.
10 mcor[5,2] -0.501 -0.515 0.137 0.134 -0.701 -0.251 1.00 3243. 3085.
# … with 15 more rows
sdat <- list(n=n,d=d, dat=dat,corpriorscale= 0.2,
x = 0.9)
> mod_out$summary("mcor")
# A tibble: 25 x 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 mcor[1,1] 1 1 0 0 1 1 NA NA NA
2 mcor[2,1] 0.902 0.907 0.0345 0.0316 0.836 0.950 1.00 1397. 1766.
3 mcor[3,1] 0.611 0.616 0.0889 0.0892 0.458 0.750 1.00 2151. 2387.
4 mcor[4,1] 0.718 0.725 0.0777 0.0756 0.575 0.833 1.00 1771. 2087.
5 mcor[5,1] 0.556 0.559 0.0937 0.0955 0.400 0.701 1.00 1964. 2642.
6 mcor[1,2] 0.902 0.907 0.0345 0.0316 0.836 0.950 1.00 1397. 1766.
7 mcor[2,2] 1 1 0 0 1 1 NA NA NA
8 mcor[3,2] 0.720 0.723 0.0673 0.0672 0.605 0.824 1.00 2942. 2419.
9 mcor[4,2] 0.942 0.946 0.0247 0.0213 0.894 0.973 1.00 2253. 2528.
10 mcor[5,2] 0.608 0.612 0.0779 0.0787 0.473 0.729 1.00 2801. 2709.
```