How to sample from a sum of two distributions: binomial and poisson

Question

Is there a way to predict a value from a sum of two distributions? I am getting a syntax error on rstan when I try to estimate y here: y ~ binomial(,) + poisson()

library(rstan)

BH_model_block <- "
data{
  int y; 
  int a; 
}

parameters{
  real <lower = 0, upper = 1> c;
  real <lower = 0, upper = 1> b;
}

model{
  y ~ binomial(a,b)+ poisson(c);
}
"
BH_model <- stan_model(model_code = BH_model_block)
BH_fit <- sampling(BH_model,
                   data = list(y = 5,
                               a = 2), 
                   iter= 1000)

Produces this error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

  error in 'model2c6022623d56_457bd7ab767c318c1db686d1edf0b8f6' at line 13, column 20
  -------------------------------------------------
    11: 
    12: model{
    13:   y ~ binomial(a,b)+ poisson(c);
                           ^
    14: }
  -------------------------------------------------

PARSER EXPECTED: ";"
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model '457bd7ab767c318c1db686d1edf0b8f6' due to the above error.

Answer 1

An alternative is to substitute the Binomial with a Poisson, and use Poisson additivity:

BH_model_block <- "
data{
  int y; 
  int a; 
}

parameters{
  real <lower = 0, upper = 1> c;
  real <lower = 0, upper = 1> b;
}

model{
  y ~ poisson(a * b + c);
}
"

This differs in that if b is not small, the Binomial has a lower variance than the Poisson, but maybe there is overdispersion anyhow?