# Sum of independent beta-Bernoulli random variables

#### SBlais

##### New Member
I'm looking for the distribution of the sum of independent beta-Bernoulli random variables, that is Bernoulli random variables with beta-distributed parameters, x_i ~ Bern(p_i) with p_i ~ beta(alpha,beta) (same alpha and beta, of course). I guess one can work with convolutions or characteristic functions.

The distribution would have three parameters: number of trials, alpha and beta.

Another way to approach the problem would be to consider a beta-mixture of Poisson-Binomial random variables and integrate the p_i's out.

I'm hoping someone somewhere has already done the math, so I don't have to do it myself...

#### BGM

##### TS Contributor
Actually a Beta-Bernoulli random variable is a Bernoulli random variable so their independent sum is just a Binomial distribution.

Not sure if this is what you want.

#### SBlais

##### New Member
Actually a Beta-Bernoulli random variable is a Bernoulli random variable so their independent sum is just a Binomial distribution.

Not sure if this is what you want.
Not exactly, but thanks . That would be true if the parameters of the Bernoulli were fixed and the same, i.e. sum of Bern(p)'s. Here, the parameters are different and random, but from the same distribution, i.e. sum of Bern(p_i) with p_i ~ beta(alpha, gamma). In other words, I have a population of experiments, each of which has its own probability of success, p_i. The population is such that the p_i's are beta(alpha,beta). I want the distribution of the number of successes for a sample of size N from that population. The distribution will have three parameters: N, alpha and beta. (A binomial distribution assumes p is known and has two parameters, N and p.)

#### Dason

Yeah - each independent trial is coming from a bernoulli-beta distribution. But a bernoulli-beta distribution is just a bernoulli because you can integrate out alpha and beta.

Have some R code if you don't believe us:
Code:
# Set alpha and beta
alp <- 5
bet <- 7

# How many of these beta-bernoullis do we want to sum?
n <- 10
# How many trials do we want to do?
N <- 100000

# Generate the success probability from the beta
pi <- rbeta(n*N, alp, bet)
# Generate the observation from the binomial
x <- rbinom(n*N, 1, pi)
# Get the desired sum
y <- apply(matrix(x, ncol = n, byrow = T), 1, sum)

# We claim that obs will follow a binomial with n = n, p = alp/(alp + bet)
# This gets us the empirical estimate of the pmf
obs <- table(y)/length(y)
# This gets us the theoretical estimate of the pmf (but only the ones we got outcomes for)
exp <- dbinom(as.numeric(names(table(y))), n, alp/(alp + bet))

# Look at them together.  If these are essentially the same that provides evidence
# that it just ends up being binomial
cbind(obs, exp)
Results:
Code:
> cbind(obs, exp)
obs          exp
0  0.00456 0.0045621330
1  0.03279 0.0325866641
2  0.10397 0.1047428489
3  0.19746 0.1995101884
4  0.25227 0.2493877355
5  0.21331 0.2137609161
6  0.12771 0.1272386405
7  0.05132 0.0519341390
8  0.01422 0.0139109301
9  0.00228 0.0022080841
10 0.00011 0.0001577203

#### SBlais

##### New Member
We claim that obs will follow a binomial with n = n, p = alp/(alp + bet)
There's an extra bit of information here: binomial with p = alpha / (alpha+beta). That's what I was looking for (it is correct, indeed). Thanks a lot.

Yeah - each independent trial is coming from a bernoulli-beta distribution. But a bernoulli-beta distribution is just a bernoulli because you can integrate out alpha and beta.

But, by the way, you're not integrating alpha and beta out here, since they still show up in the distribution... it is the p's that are integrated out. That was actually the second approach to solve the problem that I had mentionned: to intergrate out the p's from a beta-mixture of Poisson-binomial. It seems this is a classical result in your field - I'm from time series. Do you have a reference to propose?

Thanks again.