I have the stan model below implementing a multivariate normal model (the parameters are scaling variance and regression constants) with a custom log likelihood.

I am working on implementing this with a separable covariance (i.e Sigma can be factored in A \otimes B).

As part of this I need to take the residual vector Y-XB (here called mu_minus_y) and stack it into a matrix with column dimension equal to the dimension of B and row dimension equal to the dimension of A.

I am looking at the rstan documentation to see if there is a built in function I can use for this but have been unable to find one.

Does anyone happen to know a function I can use for this?

If not, will I need to build this as a c++ function or can I build a custom function using other stan built ins?

Thank you

```
functions{
real new_multiNormal_lpdf(vector y, vector mu, matrix sigma,int n){
vector[n] mu_minus_y_;
vector[n] quadratic_form_1_;
real quadratic_form_;
real log_det_;
real log_lik_;
real lprob_new_normal_;
mu_minus_y_=y-mu;
quadratic_form_1_= mdivide_left_spd(sigma,mu_minus_y_);
quadratic_form_=mu_minus_y_'*quadratic_form_1_;
//
log_det_=log_determinant(sigma);
//lprob_new_normal_=multi_normal_lpdf(y|mu,sigma);
lprob_new_normal_= -.5*log_det_-.5*quadratic_form_ -log(2*pi());
return lprob_new_normal_;
}
}
data {
int<lower = 1> n;
int<lower = 1> p;
matrix[n, p] X;
matrix[n,n] sigma;
vector[p] beta_0;
matrix[p,p] sigma_0;
//vector[n] mu;
vector[n] y;
}
parameters {
vector[p] beta;
real<lower = 0> rho;
}
model {
beta ~ multi_normal(beta_0, sigma_0);
rho ~ inv_gamma(2, 2);
//y ~ multi_normal(X*beta, rho*sigma);
target+=new_multiNormal_lpdf(y|X*beta,rho*sigma,n);
}
```