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.

  • 6
    $\begingroup$ "I was getting results of ~1000 samples needed for a maximum error of 3% at 95% confidence." --- pretty much; that's why polls often sample about 1000 people ... and then report a margin of error in the order of 3%. It applies pretty well when the percentage isn't close to 0 or 1 (it's too wide in those cases) $\endgroup$
    – Glen_b
    Jan 6 '15 at 23:47
  • 1
    $\begingroup$ What do you mean by "which even is more efficient"? Do you mean "which event has the larger expected reward"? $\endgroup$
    – Glen_b
    Jan 6 '15 at 23:50
  • 1
    $\begingroup$ Yeah, I'm trying to figure out which has the larger expected reward over time. I can't do both events -- have to choose one or the other. While Event #1 gives less reward, it is possible it gives the reward more often. $\endgroup$
    – Brad
    Jan 7 '15 at 16:05
  • 2
    $\begingroup$ You could use sequential samplig, in which sample size is not fixed in advance. The advantage to this approach is that it guarantees a confidence no matter what the (unknown) probability be. See for example here; specially the last referenced paper $\endgroup$
    – Luis Mendo
    Jan 9 '15 at 16:06

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.


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)

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)

     temp <- store[idx_p,]

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

     temp2 <- temp[idx_2,]

     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
     xlab = "prevalence", ylab = "sample size")

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.

sample size vs prevalence

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.

  • $\begingroup$ If this is slow, it is because you make the steps way too small! $\endgroup$ Nov 7 '18 at 9:35
  • $\begingroup$ @kjetilbhalvorsen - it was good enough. The "slow" was a disclaimer, and a "handrail" that can help the asker feel more comfortable engaging the method. If you don't know the analytic approach, a simulation can help you teach yourself, or help you decide if you need to ask for help. $\endgroup$ Nov 7 '18 at 14:04

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.


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.