# Sum of N binomials with Stan

I'm having trouble implementing a sum of $$N$$ binomials (and a poisson distribution) with Stan.

Observed data is $$(y_i)$$ for $$1 \leq i \leq M$$, and $$(x_{ij})$$ for $$1 \leq i \leq M$$, $$1 \leq j \leq N$$.

Assuming the following generating process :

$$y_i = \sigma + \sum_{j=1}^N y_{ij} \\ \textrm{with } y_{ij} \sim \textrm{Binomial}(x_{ij},p_j)\\ \textrm{and } \sigma \sim \textrm{Poisson}(\lambda)$$

I would like to fit my data against the paramters $$(p_j)$$ and $$\lambda$$, using informative priors on $$p_j$$ (though this is not relevant here). But I don't know how to implement this with Stan. Binomials are discrete and only produce integers which are not allowed as Stan parameters. Besides, I have no idea how to model a sum of N Binomials + a Poisson variable in Stan.

Assuming the binomials are well approximated by poisson distributions would make things much easier. We then obtain :

$$y_i \sim \textrm{Poisson}(\lambda + \sum_{j=1}^N p_{j} x_{ij})$$

I wrote the following stan-model based on this approximation :

data {
int<lower=0> M;
int<lower=0> N;
int<lower=0> y[M];
int<lower=0> x[M,N];
}
parameters {
vector <lower=0,upper=1>[N] p;
real <lower=0> sigma;
}

transformed parameters {
real<lower=0> yf[M];

for (i in 1:M) {
yf[i] = sigma;
for (k in 1:N) {
yf[i] += p[k]*x[i,k];
}
}
}

model {
p ~ uniform(0,1); // un-informative prior for simplification
sigma ~ uniform(0,50);
y ~ poisson(yf);
}


In my case, $$x$$'s orders of magnitude range between $$0$$ and $$10^5$$, and $$p$$'s between $$10^{-6}$$ and $$10^{-2}$$. Usually $$x_{ij}p_j \sim \mathcal{O}(1)$$ but sometimes it can exceed a few hundreds. Therefore, I doubt the Poisson approximation is really appropriate.

Is Stan suited for this problem ? If it is, how to implement this model without requiring the Poisson approximation ?

Stan is able to fit such a model. I guess that in your model $$\sum_{j=1}^{N} y_{i j}$$ is there to capture the mean and $$\sigma$$ to capture the deviation. I would then propose the following adaptation:
$$y_{i} \sim \mathcal{N}(\sum_{j=1}^{N} y_{i j}, \sigma)$$
$$\text { with } y_{i j} \sim \operatorname{Binomial}\left(x_{i j}, p_{j}\right)$$
$$\text { and } \sigma \sim \operatorname{Poisson}(\lambda)$$
In the model block, you can say that the mean is the sum of $$y_{ij}$$ and give a binomial prior to your $$y_{ij}$$.