Bootstrapping probability of success: what is wrong with my method? I have a dataset where each row represents the outcome of individual bets placed by a population of players; each row has two fields: a player identifier and
the outcome of the bet (1 if the player won, 0 if he lost). 
On my laptop I can not load all the dataset into memory, so I have created
a different dataset (that fits into memory) having a row for each player made of the number of bets placed by the player (different for each player) and the number of bets he won.
I can now calculate the probability of winning of the population by summing the bets won by all the players, summing the number of bets placed by the players and estimating the ratios, let's call it p.
I want to know how variable is p, so I applied the bootstrapping technique to the samples: for one thousand time I resample with replacement the players,
calculate the total number of bets won, calculate the total number of bets placed and estimate the ratios. So in the end I get 1000 values p, for which I can look at the distribution.
If I look at the distribution of p obtained with the bootstrapping I get a different distribution than the one I obtain simulating the process drawing 1 and 0 with the value of p estimated by the whole population.
I know that the second process is the correct one (I also know that for large n I the distribution of p is normal, and I can calculate mu and sigma directly from my dataset), but what I really would like to know is why bootstrapping p
with the method described above generates a distribution that differs from the theoretical one, or by the one obtained simulating the process.
 A: Still not sure if I understood everything what you are trying to do, but I will try adding some comments.
Trying to repeat your experiment
First I tried to replicate your experiment (so you can check if it's similar to what you are actually doing).
I created a table of 10 players:
playerBets <- Map(rbinom, n=runif(10, min=1, max=100), 1, runif(10))
names(playerBets) <- letters[1:10]

theData <- cbind(sapply(playerBets, length), sapply(playerBets, mean))
theData[,2] <- theData[,2]*theData[,1]

theData
   [,1] [,2]
 a   68   47
 b   10    0
 c   62   19
 d   39   13
 e   78    1
 f   41   21
 g   61   28
 h   91   37
 i    7    3
 j   79   66

Here first column in number of bets, 2nd column - number of wins.
Then a simple function to get p scores:
getP <- function(myData) {
    sum(myData[,2]) / sum(myData[,1])
}

estimated p:
p <- getP(theData)
p
[1] 0.4384328

Bootstrap by sampling the players.
Here I simply draw 10 rows from my data with replacement (since I have 10 players)
boot_ps <- replicate(1000, {
    getP(theData[sample(nrow(theData), replace=TRUE), ])
})
1000 simulations using real p:
Here I am drawing from binomial distribution (as many times as there were total bets in the data, using estimated p). And repeating this 1000 times as well to obtain a distribution.
simul_ps <- replicate(1000, mean(rbinom(sum(theData[,2]), p=p, size=1)))

Plot densities:
plot(density(simul_ps), type="l", col="green", xlim=c(0,1))
points(density(boot_ps), type="l", col="red", xlim=c(0,1))
legend("topright", legend=c("simulation","bootstrap"), c("green","red"), c("green","red"))


Comments
So, assuming this is similar to what you get...
When you bootstrap players from the table - you take all the info in one block - their p and their n.
A few things can happen in this case:

First

your bootstraps all have different n - number of trials. Very simple example: you have 2 players. One played 1000 times, and one played 2 times. p for the first is 0.25, p for the second is 1 (since he won all 2 by accident). Now when you bootstrap you will be getting some scenarios when all you draw is the second player 2 times. so you will have p = 1, but the n in this case is only 4.

Second

Your players come in blocks. So assume you have one big winner and he played a 1000 games. In reality the number of games he plays can vary. However in your simulation you can only change that by thousands (draw 1000 from such a player, draw 2000, draw 3000, etc).
Modified bootstrap
The way I would think about this experiment is something like this:


*

*Your population is composed of bets that have winning probabilities. $p(win)$

*You also have different types of players. $pl$

*So each bet has winning probability $p(win)$, which is conditioned on the type of player that made it. $p(win|player=pl)$

*Players have different probabilities of actually making a bet $p(pl)$. So then each bet has a joint probability of: what kind of player made a bet, what is the probability of that player making a bet and what is the probability of that player winning. $p(win|player=pl)p(pl)$


So when you do simulation, you can do it this way:


*

*Estimate p for all players. So you get $p_{player}$ for each player.

*Make a populations of p estimates by repeating each player's p the number of times that player made a bet. Thus now you get $p_{player}n_{player}$ where $n_{player}$ is simply the number of times that player made a bet.

*Now bootstrap those p scores with replacement. You will obtain a population of players making bets. So you will obtain a new sample of players from the population. Doing this you will be changing $n_{player}$ into $p(player)$ by via bootstrapping. You get $p_{player} * p(player)$

*Finally draw from binomial distribution, using bootstrapped p for each draw. That is - you now simulate each player making a bet again. Each time a player can loose/win with his probability p. This will give you $p(win|player)p(player)$ over all players.

*Repeat 3-4 1000 times.


If I understood your bootstrapping correctly then you did bootstrap on the values obtained in the second step. Notice how your probabilities and number of bets come in blocks. And both the probability and number of bets are fixed for one block (they are not being permuted).
And your other simulation using p tried to estimate $p(win)$ It had the same number of trials (number of bets made) in each case, but it assumed that all players had the same probability of winning.
Here is how I did the modified version:
thePs <- rep(theData[,2]/theData[,1], theData[,1])
boot2_ps <- replicate(1000, mean(rbinom(sum(theData[,1]), p=sample(thePs, replace=TRUE), size=1)))

This is what I got:

Please note that I am not saying that the modified version is the correct way to do it. But it should make the illustrations and comments more clear.
