Probability theory can be a weird place sometimes. Here I was, confident in my insane math skills, trying to solve the following problem:

Let $N, \alpha$ and $\beta$ be given.

  • Draw $N$ coins from a fixed bag with different coins, such that the success probability for each coin is unknown but modeled by a Beta distribution $Beta(\alpha, \beta)$.
  • Flip each coin from the bag once and count the number of successes.

How are the number of successes distributed?

Realisation of ignorance

"Pff, what an insult!", I thought. "It's obviously the Beta-Binomial distribution!", my inner Icarus cried. You can simply check that with an easy (Python) implementation:

sample_size = 2**20

N = 42
alpha = 3.14159
beta  = 2.71828

coins = np.random.beta(alpha, beta, size=N)

sample = []
for i in range(sample_size):
    R = np.random.uniform(size=N)
    wins = np.sum(np.where(coins>R, 1, 0))

Now we can compare mean and variance and that should be enough for the sanity check:

>>> print(np.mean(sample), N*alpha/(alpha+beta))
22.00448989868164 22.517014882582718

Well, that error is tolerable; we're dealing here with probabilities after all. Let's move on to the variance

>>> print(np.var(sample), N*alpha*beta*(alpha+beta+N)/((alpha+beta)**2 * (alpha+beta+1)))
8.996273862416274 72.87400738804808


This can't be true. Lemme draw a histogram:

enter image description here

A desperate attempt

Okay. After sitting down for a while, questioning my existence, I swallow my pride and complied with what I perceived as reality. The simulated experiment kinda looked like a binomial distribution, let's do a quick check:

enter image description here

where $\widehat{\mu}$ is the sample mean. That doesn't look too bad after all, maybe a $\chi^2$-test can formalize our visuals:

>>> cats, f_obs = zip(*[(cat, len(list(group))) for cat,group in groupby(sorted(sample))])
>>> f_exp = scipy.stats.binom.pmf(k=cats, n=N, p=np.mean(sample)/N)*sample_size
>>> scipy.stats.chisquare(f_obs=f_obs, f_exp=f_exp)
Power_divergenceResult(statistic=14303.936045316765, pvalue=0.0)

*drops the mic and leaves the room in shame*

Cry for help

At this point I gave up completely on making any further pseudo-educated guesses and started to consult the Interweb for its glorious wisdom. I stumbled upon this answer but @probabilityislogic seems to be suggesting to do exactly what I have done.

Can someone please help me to get some insights for the distribution that results from the above mentioned experiment?


1 Answer 1


Your error is that you use precisely the same $42$ probabilities (coins) in each of the $2^{20}$ simulations, and so lose the element of sample variance which would come from these varying

You would get the Beta-Binomial distribution if you chose new probabilities each time, for example putting

coins = np.random.beta(alpha, beta, size=N)

inside the for loop

  • $\begingroup$ Thanks for this insight! I guess my question then is not well posed. In my actual problem where I need this the bag of coins indeed is fixed, and we only know that the success probabilities of the (unknown but fixed) coins are modeled by a Beta-distribution. $\endgroup$ Feb 22, 2019 at 10:45
  • $\begingroup$ I modified the question box now. $\endgroup$ Feb 22, 2019 at 10:48

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.