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

user_6664233 0
Viewed 321

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.
rstan bayesian
2 Answers

Stan doesn't support integer parameters, so you can't technically do that. For two real variables, it'd look like this:

parameters {
  real x;
  real y;
}
transformed parameters {
  real z = x + y;
}
model {
  x ~ normal(0, 1);
  y ~ gamma(0.1, 2);
}

Then you get the sum distribution for z. If the variables are discrete, it won't compile.

If you don't need z in the model, then you can do this in the generated quantities block,

generated quantities {
  int x = binomial_rng(a, b);
  int y = poisson_rng(c);
  int z = x + y;
}

The drawback of doing this is that none of the variables are available in the model block. If you need discrete parameters, they need to be marginalized as described in the user's guide chapter on latent discrete parameters (also in the chapter on mixtures and HMMs). This is not so easy with a Poisson, because support isn't bounded. If the expectations of the two discrete distributions is small, then you can do it approximately with a loop over plausible values.

It looked from the example in the original post that z is data. That's a slightly different marginalization over x and y, but you only sum over the x and y such that x + y = z, so the combinatorics are greatly reduced.

user_1541470

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?

user_9332280
Related