Estimating successes while obtaining Bernoulli samples

Question

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

Answer 1

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)
[1] 0.5

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

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

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

Answer 2

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