Predict in generated quantities for mixture models

Modeling
1

I use stan to fit zero-inflated beta regression, which is mixture models of logistic regression and beta regression,

 for (n in 1:N) {
    if (y[n] == 0) {
      target += bernoulli_logit_lpmf(1 | zi);
    } else {
      target += bernoulli_logit_lpmf(0 | zi) 
                  + beta_lpdf(y[n] | mu[n] * phi[n], (1 - mu[n]) * phi[n]);
    }
  }

The complete code is as follows

data {
  int<lower=0> N;                            
  int<lower=0> K;                            
  matrix[N, K] X;                     
  vector[N] y;  
}
parameters {
  vector[K] beta;  
  vector[K] beta_phi;       
  vector[K] beta_zi;       
}
model {
  vector[N] mu  = X * beta;
  vector[N] phi = X * beta_phi;
  vector[N] zi  = X * beta_zi;

  for (n in 1:N) {
    mu[n] = inv_logit(mu[n]);
    phi[n] = exp(phi[n]);
  }
  
  for (n in 1:N) {
    if (y[n] == 0) {
      target += bernoulli_logit_lpmf(1 | zi);
    } else {
      target += bernoulli_logit_lpmf(0 | zi) 
                  + beta_lpdf(y[n] | mu[n] * phi[n], (1 - mu[n]) * phi[n]);
    }
  }
    
  target += student_t_lpdf(beta | 3, 0, 2.5);
  target += student_t_lpdf(beta_phi | 3, 0, 2.5);
  target += logistic_lpdf(beta_zi[1] | 0, 1);
  
}

The code above works fine. My question is how to write code in generated quantities { } to predict on new data, such as X_new，I do not know how to mix the bernoulli_logit_rng() and beta_rng()

generated quantities {
  vector[M] y_predict; 

  vector[M] mu  = X_new * beta;
  vector[M] phi = X_new * beta_phi;
  vector[M] zi  = X_new * beta_zi;
  
  for (i in 1:M) {
    mu[i] = inv_logit(mu[i]);
    phi[i] = exp(phi[i]);
  }
    
  // how to code ?
  for(i in 1:M) {
     //  bernoulli_logit_rng(zi);
     //  beta_rng(mu[i] * phi[i], (1 - mu[i]) * phi[i]);
    }


 }

Thanks for your help.

2

I don’t know if it’s correct

data {
  int<lower=0> N;                            
  int<lower=0> K;                            
  matrix[N, K] X;                     
  vector[N] y;  
  matrix[2, K] new_X;
}

parameters {
  vector[K] beta;  
  vector[K] beta_phi;       
  vector[K] beta_zi;       
}
model {
  vector[N] mu  = X * beta;
  vector[N] phi = X * beta_phi;
  vector[N] zi  = X * beta_zi;

  for (n in 1:N) {
    mu[n] = inv_logit(mu[n]);
    phi[n] = exp(phi[n]);
  }
  
  for (n in 1:N) {
    if (y[n] == 0) {
      target += bernoulli_logit_lpmf(1 | zi[n]);
    } else {
      target += bernoulli_logit_lpmf(0 | zi[n]) 
                  + beta_lpdf(y[n] | mu[n] * phi[n], (1 - mu[n]) * phi[n]);
    }
  }
    
  target += student_t_lpdf(beta[1] | 3, 0, 2.5);
  target += student_t_lpdf(beta_phi[1] | 3, 0, 2.5);
  target += logistic_lpdf(beta_zi[1] | 0, 1);
  
}

generated quantities {
  vector[2] y_predict; 

  vector[2] mu  = new_X * beta;
  vector[2] phi = new_X * beta_phi;
  vector[2] zi  = new_X * beta_zi;
  
  for (i in 1:2) {
    mu[i] = inv_logit(mu[i]);
    phi[i] = exp(phi[i]);
  }
  
  for(i in 1:2) {
    real tmp = bernoulli_logit_rng(zi[i]);
    if (tmp == 1) {
      y_predict[i] = 0;
    } else {
      y_predict[i] = beta_rng(mu[i] * phi[i], (1 - mu[i]) * phi[i]);
    }
 }
 
}