Constrain the mean of some variables of MVN to 0 if the MVN was reparametrized to a non-central parameterization

Modeling
1

Hello,

I am trying to fit a model which includes a multivariate normal distribution which was reparameterized to a non-central parameterization. Therefore, the sampling block doesn’t explicitly contain the multivariate normal distribution and I can’t just pass a vector of 0s for the means of the MVN variables.

This is the code which I am using:

data {
  int<lower=0> n_person;
  int<lower=0> n_item;
  int<lower = 1> Dim;
  int<lower = 1> N;
  int<lower = 1> jj[N];
  int<lower = 1> ii[N]; 
  int<lower=0, upper = 1> Y[N];
  vector<lower=0>[N] t; 
  vector<lower=0, upper=1>[N] X; 
  vector[Dim-1] Zero; 
  vector<lower=-200, upper = 200>[n_person] M;
  int<lower = 1> n[n_person];
  vector[n_item] b;
  vector[n_item] beta;
}

parameters {
  cholesky_factor_corr[Dim] choleskyP;
  matrix[n_person, Dim] PersParStar;
  vector<lower=0,upper=pi()/2>[Dim] sigmaPStar;
  real<lower=0,upper=pi()/2> sigmaTStar;
  vector[n_person] gamma0;
  vector[n_person] gamma1;
  vector[n_person] delta0;
  vector[n_person] delta1;
}

transformed parameters {
  vector<lower=0>[Dim] sigmaP = 2.5 * tan(sigmaPStar);
  real<lower=0> sigmaT = 2.5 * tan(sigmaTStar);
  matrix[n_person, Dim] PersPar = PersParStar * diag_pre_multiply(sigmaP, choleskyP)';
  matrix[Dim, Dim] correlP = multiply_lower_tri_self_transpose(choleskyP);
  real theta0 = mean(col(PersPar, 1));
  real zetaTheta1 = mean(col(PersPar, 2));
  real tau0 = mean(col(PersPar, 3));
  real zetaTau1 = mean(col(PersPar, 4));
  correlP[2,5] = 0;
  correlP[4,5] = 0;
  correlP[5,2] = 0;
  correlP[5,4] = 0;
  PersPar[n,5] = M[n];
  theta0 = 0;
  zetaTheta1 = 0;
  tau0 = 0;
  zetaTau1 = 0;
}

model {
  target += std_normal_lpdf(to_vector(PersParStar));
  choleskyP ~ lkj_corr_cholesky(1);
  b ~ normal(0, 2);
  beta ~ normal(0, 2);
  sigmaPStar ~ uniform(0, pi()/2);
  sigmaTStar ~ uniform(0, pi()/2);
  target += bernoulli_logit_lpmf(Y | (PersPar[ii,1] + (gamma0[ii] + gamma1[ii] .* PersPar[ii,5] + PersPar[ii,2]) .* X) - b[jj]);
  target += lognormal_lpdf(t | beta[jj] - (PersPar[ii,3] + (delta0[ii] + delta1[ii] .* PersPar[ii,5] + PersPar[ii,4]) .* X), sigmaT);
}

My attempt to constrain the means to 0 in the transformed parameters block is not working, since after the estimation has been done and I have the results the means of the MVN variables are -0.5, -10.93, -7.42 and 1.3. How can I constrain these means to always be 0?