// glmMixBayes_binomial.stan

data {
  // Core Data
  int<lower=1> N;               // Number of observations
  int<lower=1> K;               // Number of predictors
  matrix[N, K] X;               // Predictor matrix
  int<lower=0, upper=1> y[N];   // Binary outcome variable (0 or 1)

  // Prior Hyperparameters (Passed from R)
  // Assuming normal() for betas and beta() for theta.
  real prior_beta1_mu;              // e.g., 0
  real<lower=0> prior_beta1_sd;     // e.g., 2.5
  real prior_beta2_mu;              // e.g., 0
  real<lower=0> prior_beta2_sd;     // e.g., 5.0

  real<lower=0> prior_theta_alpha;  // e.g., 1.0
  real<lower=0> prior_theta_beta;   // e.g., 1.0
}

parameters {
  // Model Parameters
  real<lower=0, upper=1> theta;     // Mixing proportion (for component 1)
  vector[K] beta1;                  // Regression coefficients for component 1
  vector[K] beta2;                  // Regression coefficients for component 2
}

transformed parameters {
  // Vectorized Linear Predictors
  // Pre-computing the logit (eta) once per iteration
  vector[N] eta1 = X * beta1;
  vector[N] eta2 = X * beta2;
}

model {
  // Priors
  beta1 ~ normal(prior_beta1_mu, prior_beta1_sd);
  beta2 ~ normal(prior_beta2_mu, prior_beta2_sd);
  theta ~ beta(prior_theta_alpha, prior_theta_beta);

  // Likelihood
  for (n in 1:N) {
    target += log_mix(
      theta,
      bernoulli_logit_lpmf(y[n] | eta1[n]),
      bernoulli_logit_lpmf(y[n] | eta2[n])
    );
  }
}

generated quantities {
  // Posterior Predictions / Mixture Assignments
  int<lower=1, upper=2> z[N];
  for (n in 1:N) {
    vector[2] lw;
    lw[1] = log(theta) + bernoulli_logit_lpmf(y[n] | eta1[n]);
    lw[2] = log1m(theta) + bernoulli_logit_lpmf(y[n] | eta2[n]);
    z[n] = categorical_rng(softmax(lw));
  }
}
