This is the example of a multivariate model from the users guide :

```
data {
int<lower=1> K;
int<lower=1> J;
int<lower=0> N;
vector[J] x[N];
vector[K] y[N];
}
parameters {
matrix[K, J] beta;
cholesky_factor_corr[K] L_Omega;
vector<lower=0>[K] L_sigma;
}
model {
vector[K] mu[N];
matrix[K, K] L_Sigma;
for (n in 1:N)
mu[n] = beta * x[n];
L_Sigma = diag_pre_multiply(L_sigma, L_Omega);
to_vector(beta) ~ normal(0, 5);
L_Omega ~ lkj_corr_cholesky(4);
L_sigma ~ cauchy(0, 2.5);
y ~ multi_normal_cholesky(mu, L_Sigma);
}
```

My model is something like this, but for each of my K outcomes Y, I have a baseline value Yb. So my model looks something like this, where the new bits are **ed

```
data {
int<lower=1> K;
int<lower=1> J;
int<lower=0> N;
vector[J] x[N];
vector[K] y[N];
**vector[K] yb[N];**
}
parameters {
matrix[K, J] beta;
**vector[K] beta_yb**;
cholesky_factor_corr[K] L_Omega;
vector<lower=0>[K] L_sigma;
}
model {
vector[K] mu[N];
matrix[K, K] L_Sigma;
for (n in 1:N)
mu[n] = beta * x[n] + yb[n]*beta_yb; // ????
L_Sigma = diag_pre_multiply(L_sigma, L_Omega);
to_vector(beta) ~ normal(0, 5);
L_Omega ~ lkj_corr_cholesky(4);
L_sigma ~ cauchy(0, 2.5);
y ~ multi_normal_cholesky(mu, L_Sigma);
}
```

I am struggling with what goes before the “???”. I’ve tried every combination I can think of to produce a column vector that looks like yb[n]*beta_yb, but they all give errors, either when parsing or when trying to initialize any chain. I have tried declaring both as every combination of matrix and vector or row_vector that makes sense to combine, but nothing will get past the first iteration of any chain.

I also tried doing this manually,

```
vector[K] product[N];
for (n in 1:N) {
for (k in 1:K) {
product[n,k]=beta_yb[k]*yb[n,k];
}
}
for (n in 1:N)
mu[n] = beta * x[n] + product[n];
```

with similar (non) results. Does anyone have any thoughts on how to do this?