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 ?


1 Answer 1


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}$.

Hope it helps.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.