Binomial Probability: State Lottery Example

Each month, a state of 4,000,000 citizens is giving an iPhone to 1,000 citizens.

I would like to know the probability of somebody (any single person) winning exactly 2 times in 11 months and the probability of somebody winning exactly 3 times in 11 months.

I've seen the binomial pmf, but I can't figure out where to input the 1,000. Does it go in $p$ (probability of success in one trial) or multiply the result $P$ by 1,000?

• Hint: what proportion of people win? Protip: Use the normal approximation since the numbers are large. Commented Jan 12, 2018 at 19:32
• Please clarify the event of interest: would it be the probability that you win either 2 or 3 times in 11 months (this chance is just a few per million), or is it the probability that somebody wins 2 or 3 times in 11 months (a chance of about 13%)?
– whuber
Commented Jan 12, 2018 at 20:20
• I would like to know the latter Commented Jan 16, 2018 at 10:40

I did some simple simulations (in R) which can give some indications. But simulating this in R takes quite some time, so to get more replications we need a more efficient implementation ... Following the code:

# For this algorithm,  see B Ripley Stochastic Simulation page 80
reservoir_sampling <- function(N, pods) {
stopifnot(N > pods)
csample <- 1:pods
for (seen in seq(from=pods+1, to=N, by=1))if(runif(1)<=(pods/seen))csample[sample(1:pods, 1)] <- seen
return(csample)
}

sim_one <- function(N=4000000, k=11, pods=1000) {
winners <- matrix(0L, nrow=k, ncol=pods)
for (round in seq(from=1, to=k, by=1)) winners[round, ]<-reservoir_sampling(N,  pods)
wintab <- table(c(winners))
names(wintab) <- NULL
wintab <- factor(wintab, levels=1:k)
wintab <- table(wintab)
wintab
}

SIMS <- replicate(1000, sim_one() )


In this 1000 replications I get somebody winning twice in ALL the replicas, somebody winning thrice in 11 replicas, so quite far from @whubers value in comment of about 13%.