I am trying to estimate a random intercept zero inflated poisson model. The non centered parametrization seems to work fine as in completes estimation but is not able to retrevie the true parameters, however, the centered parametrization does not mix even at adapt_delta set to 0.99 and 10000 iterations. It also throws low BFMI errors

Here is the code I am using for the centered parametrization:

```
data {
int<lower=0> N;
int<lower=0> K; // Number of predictors in zero inflation model - homogenous
int<lower=0> J;// Number of predictors in count model - homogenous
matrix[N,K] x1;// data for logit model - homogenous
matrix[N,J] x2;// data for count model - homogenous
int<lower=0> y[N];//count variable
int<lower=1> L;//Number of customers/groups
int<lower=1,upper=L> ll[N];//Number of observations in groups
}
transformed data {
matrix[N, K] Q_ast1;
matrix[K, K] R_ast1;
matrix[K, K] R_ast_inverse1;
matrix[N, J] Q_ast2;
matrix[J, J] R_ast2;
matrix[J, J] R_ast_inverse2;
// thin and scale the QR decomposition
Q_ast1 = qr_thin_Q(x1) * sqrt(N - 1);
R_ast1 = qr_thin_R(x1) / sqrt(N - 1);
R_ast_inverse1 = inverse(R_ast1);
Q_ast2 = qr_thin_Q(x2) * sqrt(N - 1);
R_ast2 = qr_thin_R(x2) / sqrt(N - 1);
R_ast_inverse2 = inverse(R_ast2);
}
parameters {
vector[K] zeta;//count model estimates - homogenous
vector[L] rbeta; //count model estimates - heterogenous
real rbeta_mu;//coefficient of variation of random effects for zero inflation model
real<lower=0> rbeta_sigma;//std dev of distribution of rbeta parameters
vector[J] omega;//zero inflation model estimates - homogenous
vector[L] rgamma;//zero inflation model estimates - heterogenous
real rgamma_mu;//mean of distribution of rgamma parameters
real<lower=0> rgamma_sigma;//std dev of distribution of rgamma
}
transformed parameters {
vector[N] theta;
vector[N] lambda;
vector[N] eta1;
vector[N] eta2;
eta1 = Q_ast1 * zeta;
eta2 = Q_ast2 * omega;
for(n in 1:N){
theta[n] = inv_logit(eta1[n] + rbeta[ll[n]]);
lambda[n] = exp(eta2[n] + rgamma[ll[n]]);
}
}
model {
zeta ~ normal(0,100);
rbeta_mu ~ normal(0,100);
rbeta_sigma ~ cauchy(0,10);
omega ~ normal(0,100);
rgamma_mu ~ normal(0,100);
rgamma_sigma ~ cauchy(0,10);
rbeta~normal(rbeta_mu,rbeta_sigma);
rgamma~normal(rgamma_mu,rgamma_sigma);
for (n in 1:N) {
if (y[n] == 0)
target += log_sum_exp(bernoulli_lpmf(1 | theta[n]),
bernoulli_lpmf(0 | theta[n])
+ poisson_lpmf(y[n] | lambda[n]));
else
target += bernoulli_lpmf(0 | theta[n])
+ poisson_lpmf(y[n] | lambda[n]);
}
}
generated quantities {
vector[K] beta;
vector[J] gamma;
beta = R_ast_inverse1 * zeta; // coefficients on x1
gamma= R_ast_inverse2 * omega;// coefficients on x2
}
```

This is the code for non centered parametrization:

```
data {
int<lower=0> N;
int<lower=0> K; // Number of predictors in zero inflation model - homogenous
int<lower=0> J;// Number of predictors in count model - homogenous
matrix[N,K] x1;// data for logit model - homogenous
matrix[N,J] x2;// data for count model - homogenous
int<lower=0> y[N];//count variable
int<lower=1> L;//Number of customers/groups
int<lower=1,upper=L> ll[N];//Number of observations in groups
}
transformed data {
matrix[N, K] Q_ast1;
matrix[K, K] R_ast1;
matrix[K, K] R_ast_inverse1;
matrix[N, J] Q_ast2;
matrix[J, J] R_ast2;
matrix[J, J] R_ast_inverse2;
// thin and scale the QR decomposition
Q_ast1 = qr_thin_Q(x1) * sqrt(N - 1);
R_ast1 = qr_thin_R(x1) / sqrt(N - 1);
R_ast_inverse1 = inverse(R_ast1);
Q_ast2 = qr_thin_Q(x2) * sqrt(N - 1);
R_ast2 = qr_thin_R(x2) / sqrt(N - 1);
R_ast_inverse2 = inverse(R_ast2);
}
parameters {
vector[K] zeta;//count model estimates - homogenous
vector[L] beta_raw; //count model estimates - heterogenous
real rbeta_mu;//coefficient of variation of random effects for zero inflation model
real<lower=0> rbeta_sigma;//std dev of distribution of rbeta parameters
vector[J] omega;//zero inflation model estimates - homogenous
vector[L] gamma_raw;//zero inflation model estimates - heterogenous
real rgamma_mu;//mean of distribution of rgamma parameters
real<lower=0> rgamma_sigma;//std dev of distribution of rgamma
}
transformed parameters {
vector[N] theta;
vector[N] lambda;
vector[L] rbeta;
vector[L] rgamma;
vector[N] eta1;
vector[N] eta2;
rbeta = rbeta_mu + beta_raw * rbeta_sigma; // coefficients on x1
rgamma = rgamma_mu + gamma_raw * rgamma_sigma ; // coefficients on x2
eta1 = Q_ast1 * zeta;
eta2 = Q_ast2 * omega;
for(n in 1:N){
theta[n] = inv_logit(eta1[n] + rbeta[ll[n]]);
lambda[n] = exp(eta2[n] + rgamma[ll[n]]);
}
}
model {
zeta ~ normal(0,100);
rbeta_mu ~ normal(0,100);
rbeta_sigma ~ cauchy(0,10);
omega ~ normal(0,100);
rgamma_mu ~ normal(0,100);
rgamma_sigma ~ cauchy(0,10);
beta_raw ~ normal(0,1);
gamma_raw ~ normal(0,1);
for (n in 1:N) {
if (y[n] == 0)
target += log_sum_exp(bernoulli_lpmf(1 | theta[n]),
bernoulli_lpmf(0 | theta[n])
+ poisson_lpmf(y[n] | lambda[n]));
else
target += bernoulli_lpmf(0 | theta[n])
+ poisson_lpmf(y[n] | lambda[n]);
}
}
generated quantities {
vector[K] beta;
vector[J] gamma;
beta = R_ast_inverse1 * zeta; // coefficients on x1
gamma= R_ast_inverse2 * omega;// coefficients on x2
}
```

Any suggestions as to what is going wrong? All help is appreciated.