Sample size needed to estimate probability of "success" in Bernoulli trial

Question

Suppose a game offers an event which upon completion, either gives a reward, or gives nothing. The exact mechanism for determining whether the reward is given is unknown, but I assume a random number generator is used, and if the result is greater than some hard-coded value, you get the reward.

If I want to basically reverse-engineer what value the programmers used to determine how often the reward is given (estimated 15-30%), how do I calculate the number of samples I need?

I started with the "Estimator of true probability" section here: Checking_whether_a_coin_is_fair, but I'm not certain I'm heading down the right path. I was getting results of ~1000 samples needed for a maximum error of 3% at 95% confidence.

Ultimately, here's what I'm trying to solve:

• Event #1 gives reward 1.0R, X% of the time
• Event #2 gives reward 1.4R, Y% of the time

I'd like to estimate X & Y accurately enough to determine which event is more efficient. Large sample sizes are a problem since I can only get 1 sample every 20 minutes, at most.

Answer 1

Assuming your individual trials are independent, you observe a binomial variate $$ X \sim \text{Bin}(n,p) $$ where you decide on $n$ and want to estimate $p$. Now the maximum likelihood estimator of $p$, the sample fraction $\hat{p}=X/n$ have variance $\frac{p\cdot(1-p)}{n}\le \frac1{4n}$ which is achieved for $p=\frac12$. So the standard error is $\le 1/\sqrt{4 n} = \frac1{2\sqrt{n}}$. A large-sample aprroximate confidence interval has half-width about 2 standard errors, so to keep that at maximum $0.03$, say, you have to solve $$ \frac2{2\sqrt{n}} \le 0.03 $$ which gives $n \ge 1112$. Now you can solve for other requirements for halfwidth, the same way. If you know (or is willing to assume) that $p$ is bounded away from 0.5, you can do with somewhat less observations.

Answer 2

I know it is less elegant, but I had to simulate it. Not only did I build a pretty simple simulation, but it is inelegant and slow to run. It is good enough, though. One advantage is that, as long as some of the basics are right, it is going to tell me when the elegant approach falls down.

The sample size is going to vary as a function of the hard-coded value.

So here is the code:

    #main code
    #want 95% CI to be no more than 3% from 
    # prevalence
    #expect prevalence around 15% to 30%
    #think sample size is ~1000
my_prev <- seq(from=0.15, to=0.30, 
               by = 0.002)

samp_sizes <- seq(from=400, to=800, by = 1)
samp_sizes

N_loops <- 2000

store <- matrix(0,  nrow = 
    length(my_prev)*length(samp_sizes),
                ncol = 3)
count <- 1

#for each prevalence
for (i in 1:length(my_prev)) {

     #for each sample size
     for(j in 1:length(samp_sizes)){

          temp <- 0

          for(k in 1:N_loops){

               #draw samples
               y <- rbinom(n = 
                     samp_sizes[j],
                           size = 1,
                           prob = 
                            my_prev[i])

               #compute prevalence, store
               temp[k] <- mean(y)

          }

          #compute 5% and 95% of temp
          width <-  diff(quantile(x = temp, 
             probs = c(0.05,0.95)))

          #store samp_size, prevalence, and 
          # CI half-width
          store[count, 1] <- my_prev[i]
          store[count, 2] <- samp_sizes[j]
          store[count, 3] <- width[[1]]

          count <- count+1
     }

}


store2 <- numeric(length(my_prev))

#go through store
for(i in 1:length(my_prev)){
     #for each prevalence
     #find first CI half-width below 3%
     #store samp_size

     idx_p <- which(store[, 1] == 
               my_prev[i], arr.ind = T)
     idx_p

     temp <- store[idx_p, ]
     temp

     idx_2 <- which(temp[, 3] <= 0.03*2, 
       arr.ind = T)
     idx_2

     temp2 <- temp[idx_2, ]
     temp2

     if (length(temp2[,3])>1){
     idx_3 <- which(temp2[, 3]==max(temp2[, 
              3]), arr.ind = T)
     store2[i] <- temp2[idx_3[1], 2]
     } else {
          store2[i] <- temp2[2]
     }


}


#plot it
plot(x=my_prev, y=store2, 
     xlab = "prevalence", 
     ylab = "sample size")
lines(smooth.spline(x=my_prev,y=store2), 
        col="Red")
grid()


And here is the plot of sample size vs. prevalence such that uncertainty in 95% CI for prevalence is as close as possible to $\pm$3% without going over it.

Away from 50%, "somewhat less observations" seem to be required, as kjetil suggested.
I think that you can get a decent estimate of prevalence before 400 samples, and adjust your sampling strategy as you go.  I don't think there should be a jog in the middle, and so you might bump N_loops up to 10e3, and bump the "by" in "my_prev" down to 0.001.

Answer 3

It seems like you want to estimate for Event #1 the value of $X$ and for Event #2 the value of $Y$. You can easily use Hoeffding's inequality to determine bounds here, or if you want additive, rather than multiplicative bounds, you can use Chernoff's bound.