Checking Stan code for zero-inflated LDA

Modeling
,
1

I have been trying to fit a zero-inflated Latent Dirichlet Allocation model in Stan, which is structured as follows: Suppose there are D documents, K topics, and V words, then for the d^{th} document, \theta^{(d)} \sim \text{Dir}(\alpha). The topic assignments are done by sampling z_{dn} \sim \text{Multi}(\theta^{(d)}). The word probabilities under each topic are sampled from a zero-inflated Generalized Dirichlet distribution (ZIGD) (see Generalized Dirichlet) as follows: \beta^{(z_{dn})} \sim \begin{cases} 0 & \ \text{with probability} \ \pi \\ \text{GD}(a,b) & \ \text{with probability} \ 1-\pi \end{cases}
where \text{GD} is a Generalized Dirichlet Distribution.
Finally, the words are generated using w_{dn}|z_{dn},\beta^{(z_{dn})} \sim \text{Multi}(\beta^{(z_{dn})}). Here, I have used the stick-breaking construction to sample from a ZIGD distribution which means that if Z_j \sim \text{ZIB}(\pi_j,a_j,b_j) for j=1,2,...,k where \text{ZIB} is zero-inflated beta distribution, then take P_1=Z_1 and P_j=Z_j \prod_{i=1}^{j-1}(1-Z_i) for j=2,3,...,k. This construction gives (P_1,...,P_k) \sim \text{ZIGD}(\pi,a,b).

I am a bit skeptical about whether my code for modeling this is correct. The posterior distributions of the parameters \beta,\theta after fitting this model do not seem to match the original values of \beta,\theta in certain simulation settings and hence my question. I am not sure where exactly am I going wrong in coding this! I have written down the Stan code as follows. Thanks for reading.

data {
  int<lower=1> K; // num topics
  int<lower=1> V; // num words
  int<lower=0> D; // num docs
  int<lower=0> n[D, V]; // word counts for each doc

  // hyperparameters
  vector<lower=0>[K] alpha;
  vector<lower=0>[V] gamma1; // this is the parameter vector **a** for GD(**a**,**b**)
  vector<lower=0>[V] gamma2; // this is the parameter vector **b** for GD(**a**,**b**)
  vector<lower=0, upper=1>[V] delta; // this is the zero-inflation parameter vector $/pi$
}
parameters {
  simplex[K] theta[D]; // topic mixtures
  vector<lower=0,upper=1>[V] zeta[K]; // zero-inflated betas
}


transformed parameters {
  vector<lower=0>[V] beta[K];
  for (k in 1:K) {
	beta[k,1] =  zeta[k,1];
  for (m in 2:V) {
    beta[k,m] = zeta[k,m]*prod(1 - zeta[k,1:(m - 1)]);  // stick breaking
  }
  }
  for (k in 1:K) {
      beta[k]=beta[k]/sum(beta[k,1:V]);  // GD construction
  }
}


model {
  for (d in 1:D) {
    theta[d] ~ dirichlet(alpha);  
  }

  for (k in 1:K) {
    for (m in 1:V) {
      if (zeta[k,m]==0){  // Zero-inflated beta likelihood
        target += bernoulli_lpmf(1 | delta[m]);
      }else{
        target += bernoulli_lpmf(0 | delta[m]) + beta_lpdf(zeta[k,m] | gamma1[m], gamma2[m]);
      }
		}
  }

  for (d in 1:D) {
    vector[V] eta;
    eta = beta[1] * theta[d, 1];
    for (k in 2:K) {
      eta = eta + beta[k] * theta[d, k];
    }
    eta = eta/sum(eta[1:V]);
    n[d] ~ multinomial(eta);  // generation of each sample
  }
}

****