# Estimating successes while obtaining Bernoulli samples

I have a process which, after fixing the values of some parameters, generates samples from a Bernoulli distribution with unknown $p$.

The value of $p$ is typically small, and what I want to do is to discover suitable values for my parameters so that $p$ is at least $0.1$. An added problem is that generating each sample (i.e. running my “experiment” once) takes some considerable amount of time.

One thing I can do is, say, fix the parameters, generate 100 samples, count the number $k$ of successes and, if $k/100 < 0.1$, try again with different parameters until I find ones that yield $k/100 > 0.1$.

However, as generating each sample takes some time, intuitively I would like to stop generating samples for a fixed set of parameters if it “doesn't look promising”. For example, say, if i've already seen 30 samples and not a single success then, according to this question, with 95% confidence the value of $p < 0.1$; so it would be reasonable to stop here and not generate the remaining 70 cases. I would like to generalise this idea.

I guess the question I really want to ask is the following:

Given that I've already seen $k$ successes on $n$ samples of a Bernoulli distribution with unknown parameter $p$, what is the probability that, if I keep sampling from the same distribution, after $N$ observations (say $N = 100$), I'll see at least $K$ successes (say $K = 10$).

Update: In case anybody is interested, these are some example values I computed with a small script in python using the answers provided below. Note that both methods are trying to compute two slightly different probabilities, you might want to read the exact details below.

      Rasmus' method                | Zen's method
k=0  k=1  k=2  k=3  k=4  k=5  | k=0  k=1  k=2  k=3  k=4  k=5
n= 0  0.90                          | 0.90
n= 1  0.81 0.99                     | 0.81 0.99
n= 2  0.72 0.97 1.00                | 0.73 0.97 1.00
n= 3  0.66 0.95 1.00 1.00           | 0.66 0.95 1.00 1.00
n= 4  0.59 0.92 0.99 1.00 1.00      | 0.60 0.92 0.99 1.00 1.00
n= 5  0.53 0.89 0.98 1.00 1.00 1.00 | 0.54 0.88 0.98 1.00 1.00 1.00
n= 6  0.48 0.85 0.97 1.00 1.00 1.00 | 0.49 0.85 0.97 1.00 1.00 1.00
n= 7  0.43 0.81 0.96 0.99 1.00 1.00 | 0.45 0.81 0.96 0.99 1.00 1.00
n= 8  0.39 0.78 0.95 0.99 1.00 1.00 | 0.41 0.78 0.94 0.99 1.00 1.00
n= 9  0.35 0.74 0.93 0.99 1.00 1.00 | 0.37 0.74 0.92 0.98 1.00 1.00
n=10  0.32 0.70 0.91 0.98 1.00 1.00 | 0.34 0.70 0.90 0.98 1.00 1.00
n=11  0.28 0.66 0.89 0.97 1.00 1.00 | 0.31 0.67 0.88 0.97 0.99 1.00
n=12  0.25 0.62 0.87 0.97 0.99 1.00 | 0.28 0.63 0.86 0.96 0.99 1.00
n=13  0.23 0.58 0.84 0.96 0.99 1.00 | 0.25 0.60 0.83 0.95 0.99 1.00
n=14  0.21 0.55 0.82 0.94 0.99 1.00 | 0.23 0.56 0.81 0.93 0.98 1.00
n=15  0.19 0.51 0.79 0.93 0.98 1.00 | 0.21 0.53 0.78 0.92 0.98 0.99
n=16  0.17 0.48 0.76 0.92 0.98 1.00 | 0.19 0.50 0.76 0.90 0.97 0.99
n=17  0.15 0.45 0.73 0.90 0.97 0.99 | 0.18 0.47 0.73 0.89 0.96 0.99
n=18  0.14 0.42 0.70 0.88 0.96 0.99 | 0.16 0.45 0.71 0.87 0.95 0.98
n=19  0.12 0.39 0.68 0.87 0.96 0.99 | 0.15 0.42 0.68 0.85 0.94 0.98
n=20  0.11 0.37 0.65 0.85 0.95 0.99 | 0.14 0.40 0.65 0.83 0.93 0.98
n=21  0.10 0.34 0.62 0.83 0.94 0.98 | 0.13 0.37 0.63 0.82 0.92 0.97
n=22  0.09 0.32 0.59 0.81 0.93 0.98 | 0.12 0.35 0.60 0.80 0.91 0.96
n=23  0.08 0.29 0.56 0.79 0.91 0.97 | 0.11 0.33 0.58 0.78 0.90 0.96
n=24  0.07 0.27 0.54 0.76 0.90 0.97 | 0.10 0.31 0.56 0.76 0.88 0.95
n=25  0.07 0.25 0.51 0.74 0.89 0.96 | 0.09 0.29 0.53 0.73 0.87 0.94
n=26  0.06 0.23 0.48 0.72 0.87 0.95 | 0.08 0.27 0.51 0.71 0.85 0.93
n=27  0.05 0.22 0.46 0.69 0.86 0.95 | 0.08 0.26 0.49 0.69 0.84 0.92
n=28  0.05 0.20 0.43 0.67 0.84 0.94 | 0.07 0.24 0.47 0.67 0.82 0.91
n=29  0.04 0.18 0.41 0.65 0.82 0.93 | 0.06 0.23 0.45 0.65 0.81 0.90
n=30  0.04 0.17 0.39 0.62 0.81 0.92 | 0.06 0.21 0.43 0.63 0.79 0.89

• You appear to describe exactly the problem inverse sampling is designed to solve. – whuber Dec 3 '13 at 16:00
• Regarding inverse sampling: this question may help you – Luis Mendo Dec 4 '13 at 0:07
• Thanks for the comments. From what I could gather reading your links (and googling a bit around), using inverse sampling I would fix the target number of successes (I guess that's the $K$ w.r.t. the names I used above), and then run the experiments until $K$ successes are seen. But this does not seem to address my problem. If $p$ is very low, the time until I see those $K$ successes is going to be very large (and maybe require more than my usual 100 samples!). What I want is to be able to decide when to stop sooner because, given what I've seen, there is little hope of having $p > 0.1$. – Juan A. Navarro Dec 4 '13 at 13:55

Given that I've already seen k successes on n samples of a Bernoulli distribution with unknown parameter p, what is the probability that, if I keep sampling from the same distribution, after N observations (say N=100), I'll see at least K successes (say K=10).

A Bayesian answer to this question is: $$\Theta\sim\mathrm{U}[0,1], \qquad X\mid\Theta\sim\mathrm{Bin}(n,\Theta), \qquad Y\mid\Theta\sim\mathrm{Bin}(N,\Theta) \, .$$ $$\begin{eqnarray} P(Y\geq K&\mid& X=k) = \sum_{m=K}^N P(Y=m\mid X=k) \\ &=& \sum_{m=K}^N \int_0^1 f_{Y\mid\Theta}(m\mid\theta) f_{\Theta\mid X}(\theta\mid k)\,d\theta \qquad\qquad\qquad\qquad\qquad (*) \\ &=& \frac{\Gamma(n+2)}{\Gamma(k+1)\Gamma(n-k+1)} \sum_{m=K}^N {N \choose m} \int_0^1 \theta^{m+k}(1-\theta)^{N+n-m-k} \, d\theta \\ &=& \frac{(n+1)!}{k!(n-k)!} \sum_{m=K}^N \left( {N \choose m} \frac{\Gamma(m+k+1)\Gamma(N+n-m-k+1)}{\Gamma(N+n+2)} \right) \\ &=& \frac{(n+1)!}{k!(n-k)!(N+n+1)!} \sum_{m=K}^N \left( {N \choose m} (m+k)!(N+n-m-k)! \right)\, . \end{eqnarray}$$ The $(*)$ equality follows from the theorem of total probability, the product rule, and the conditional independence of $Y$ and $X$, given $\Theta$: $$\begin{eqnarray} f_{Y\mid X}(m\mid k) &=& \int f_{Y,\Theta\mid X}(m,\theta\mid k) \,d\theta \\ &=& \int f_{Y\mid\Theta, X}(m\mid\theta,k) f_{\Theta\mid X}(\theta\mid k) \,d\theta \\ &=& \int f_{Y\mid\Theta}(m\mid\theta) f_{\Theta\mid X}(\theta\mid k) \,d\theta \, . \end{eqnarray}$$

• Thanks a lot for this. I'm still struggling quite a bit to understand all the steps in the derivation (I have only very basic knowledge of probability), but if I just blindly plug numbers into the formula, the computed probabilities are in line with what I would expect to see (and are also quite close to those generated using Rasmus' answer). So, for the moment I'll just use the formula as is, but I would appreciate pointers to additional reading material that would help me better understand it. – Juan A. Navarro Dec 10 '13 at 3:59
• Hi, Juan. Take a look at this book by Peter Lee amazon.com/Bayesian-Statistics-Introduction-Peter-Lee/dp/… – Zen Dec 10 '13 at 14:32
• If you have any doubts, please ask us. – Zen Dec 10 '13 at 14:36

This sounds like the perfect job for Bayesian parameter estimation. So the model you have is:

$$k \sim \text{Binom}(p,n)$$

And what you want to know is the probability of $p > 0.1$. Using Bayesian parameter estimation there is nothing stopping you from estimating $\text{prob}(p > 0.1)$ at $n=1, n=2, \dots$ and stopping when the probability $\text{prob}(p > 0.1)$ is too small ("these parameters are probably not going to work") or large enough ("These parameters are probably going to work"). The two things you have to decide is: (1) What are the probability cut-offs when you are going to stop and either keep the parameters because they seem good or toss the parameters and try some new ones. (2) A prior probability for what values of $p$ are likely before you start your experiment. For (2) a reasonable starting point could be to assume a "flat prior", that is, assume that all values of $p$ are equally likely before having seen any data. If, as you state, the typical value of $p$ seems small there might be priors that better reflect the information you have.

A great introduction to how to do Bayesian inference on binomial proportions can be found here.

I don't know if you use R but a quick function that calculates $p > 0.1$ given $k$ and $n$ would be:

bin_prob <- function(k, n) {
mean(rbeta(99999, 1 + k, 1 + (n - k)) > 0.1)
}


This however assumes that before seeing any data all possible values of $p$ are equally probable. That is, before seeing any data the probability of $p > 0.1$ is 0.9. A function that perhaps would be better calibrated to your prior information would be:

bin_prob <- function(k, n) {
mean(rbeta(99999, 1 + k, 6.6 + (n - k)) > 0.1)
}


This function assumes that prior to any data the probability of $p > 0.1$ is 50%. Here follows some sample output:

> bin_prob(k=0, n=0)
 0.5

>  bin_prob(k=0, n=1)
 0.448
>  bin_prob(k=1, n=1)
 0.829

>  bin_prob(k=0, n=4)
 0.328
>  bin_prob(k=4, n=4)
 0.998

> bin_prob(k=8, n=25)
 0.997
> bin_prob(k=0, n=25)
 0.037

• thanks a lot for your solution, I learned a lot from it and from the linked introduction to Bayesian inference. Unfortunately this solution is a bit slower to compute, as it involves simulating the beta distribution a number of times in order to estimate the desired probability. In comparison Zen's exact answer can be easily computed without having to simulate random variables. – Juan A. Navarro Dec 10 '13 at 3:50
• For sure, the function can be made much faster (it almost takes as much as 30 ms on my Dell Vostro) but I just wanted go give a piece of working code that was easy to use and play around with :) – Rasmus Bååth Dec 10 '13 at 9:24
• So, as I was reading it you really want to know the probability of p > 0.1, which my answer gives. This is also what you ask after in the beginning of your question. In the end of your question you change the question to "what is the probability of seeing 10 successes out of 100", my answer does not give this. It seems like in your case it would be more usefull to know the probability of p > 0.1, unless it is the case that you actually plan t0 run exactly a hundred trials at a later stage and it it important that you then get more than 10 successes. – Rasmus Bååth Dec 10 '13 at 9:33