Sampling a discrete parameter from a bernoulli prior

Please share your Stan program and accompanying data if possible.

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

data {
  int<lower=0> N; //number of observations
  int T; //number of timepoint
  int<lower=0> Ngroups; //number of unique subjects
  int subjects[N]; //60
  vector[N] treatment;
  int familyRichness[N];
  real IgA[N];
  matrix[N, T] X_t; //intercept per time point
  matrix[N, T] X_trt; // treatment effect per time point
  matrix[N, T] X_fam; //family richness per timepoint

// The parameters accepted by the model. Our model
parameters {
  vector[T] mux; //intercept per time point for poisson model
  vector[T] muy; //intercept per time point for normal model
  vector[T] beta; // treatment effect per time point for normal
  vector[T] alpha; // treatment effect for poisson model
  vector[2] phi;
  vector[T] delta;
  real<lower=0> p4;
  real<lower=0> sigma;
  real<lower=0> sigma_iga; // random intercept variance for iga
  real<lower=0> sigma_rich; // random intercept variance for richness
  cholesky_factor_corr[2] L_iga_rich;
  matrix[2, Ngroups] z_iga_rich; //random effect from N(0, 1) to post multiply

transformed parameters {
  real<lower=0> lambda[N];
  vector[T] phis;
  real mu[N];
  matrix[2, Ngroups] ran_effs;
  vector[2] L_sigma;
  vector[T] gamma;
  L_sigma[1] = sigma_iga;
  L_sigma[2] = sigma_rich;
  phis[1] = 0;
  phis[2] = 0;
  phis[3] = phi[1];
  phis[4] = phi[2];
  for(i in 1:4)
    gamma[i] = phi[i]*delta[i];
  ran_effs = diag_pre_multiply(L_sigma, L_iga_rich) * z_iga_rich;
  for(i in 1:N) {
      lambda[i] = exp(X_t[i, ] * mux + X_trt[i, ] * alpha + ran_effs[2, subjects[i]]);//poisson
      mu[i] =  X_t[i, ] * muy + X_trt[i] * beta + X_fam[i, ]* gamma + ran_effs[1, subjects[i]]; //normal

model {
  int p1;
  int p2;
  L_iga_rich ~ lkj_corr_cholesky(2.0);
  to_vector(z_iga_rich) ~ normal(0, 1);
  sigma ~ cauchy(0, 5);
  sigma_iga ~ cauchy(0, 5);
  sigma_rich ~ cauchy(0, 5);
  mux ~ cauchy(0, 10);
  muy ~ cauchy(0, 10);
  alpha ~ normal(0, 1);
  beta ~ normal(0, 1);
  delta ~ normal(0, 1);
  phi[1] ~ bernoulli(0.5);
  phi[2] ~ bernoulli(p4);
  p4 ~ uniform(0,1);
  for(i in 1:N) {
      familyRichness[i] ~ poisson(lambda[i]);
      IgA[i] ~ normal(mu[i], sigma);

generated quantities {

  real rho;
  matrix[2, 2] Sigma;
  //compute sigma
  Sigma = diag_pre_multiply(L_sigma, L_iga_rich * L_iga_rich');
  //compute rho
  rho = Sigma[1, 2] / (pow(Sigma[1, 1], 0.5) * pow(Sigma[2, 2], 0.5));

X_{ij} \sim dpois(\lambda_{ij}), log(\lambda_{ij}) = \mu_X[j] + b1_i + \alpha[j] Z_i and
Y_{ij} \sim dnorm(\mu_{ij},\sigma), \mu_{ij} = \mu_Y[j] + b2_i + \beta[j] Z_i + \phi[j]*\delta[j] X_{ij}.
The \phi[j]'s are the challenge here. I fixed its first two components to be 0, while the last two components are to be simulated from the Bernoulli distribution. I understand that stan does not allow for declaration of discrete parameters and one has to marginalize over discrete parameters, but how do I do this in my case? Thanks in anticipation.

Hi Olajumoke! Welcome to the Stan forum! :)

Sorry it took us so long to respond to your question.The model that you posted is not so trivial it seems. Maybe you could give a bit more context, for example about what those \phi 's are doing. Also, it would be great if you could post (a subset) of your data or data that could look like your real data for us to be able to run your model.

I think then it will become a bit easier to talk about marginalizing those Bernoulli r.v.