# Bootstrap size and probability of drawing distinct observations

Suppose we draw a sample of $k$ out of $n$ numbered balls without replacement and get a sample of $N_1$ (unique) balls. We then replace all the balls in the urn and repeat the drawing a second time to get again $N_2$ (unique) balls.

We are interested in the set of balls that end up being drawn exactly once (over the course of the two experiments). It's clear that the probability, for one of the balls, of the event "being drawn exactly once" is

$$2k/n(1-k/n)$$

which is maximized at $k^*=n/2$.

Now consider the same experiment but with replacement and we look again for the value of $k$ (denoted now $\bar{k}^*$) that maximizes the frequency of the event "being drawn exactly once". Now, the expected number of distinct balls in each draw is $1-e^{-1}\approx0.632n$ and so I would have expected $\bar{k}^*=k^*/(1-e^{-1})\approx0.79n$. Running some simulations I get that the actual location of $\bar{k}^*$ seems to be much lower, probably between $0.65n$ and $0.7n$. I have a hard time understanding where is my mistake.

More generally (second question) I'm wondering what is the distribution function of the proportion of balls that are drawn exactly once in the second experiment (the one with replacements) as a function of $k$? Through numerical experimentations, I get the following curve ($n=100$): # EDIT

here is the R code to reproduce the example above:

fx02<-function(ll,n,k){
a1<-matrix(0,n,2)
a1[sample(1:n,k,replace=TRUE),1]<-1
a1[sample(1:n,k,replace=TRUE),2]<-1
sum(rowSums(a1)==1)/n
}

ss<-(1:60)*5   #The grid of values of k for which we'll compute the probability.
a4<-matrix(NA,length(ss),2)
for(i in 1:length(ss)){
a2<-ss[i]
a3<-c(lapply(1:1000,fx02,n=100,k=a2),recursive=TRUE);
a4[i,]<-c(a2,mean(a3))
}
plot(a4,xlab="k",ylab="frequency of distinct draw")

• Are you trying to maximize the chance of an individual ball being drawn exactly once or the number of balls being drawn exactly once? Jun 22, 2013 at 0:34
• Important point. The chance. How can I re-phrase my question to make this clear? Jun 22, 2013 at 0:36
• In experiment #2, what is the chance a given ball is drawn exactly once? I am thinking it is $(2k)(1/n)(1-[1/n])^{2k-1}.$ Isn't this what you want to maximize? Jun 22, 2013 at 1:28
• @soakley -- Yes, I agree. The probability of a ball being drawn once in the two samples in Experiment 2 is identical to the probability of its being drawn once in a single sample of twice the size. Less sure of this, but I believe the maximizing value of $k$ will again be $n/2$. Jun 22, 2013 at 1:36
• @soakley: the function you propose is maximized @k=n/2 (with optimal value of $\approx.0.368$)...clearly we can do better than having $.37$ chance of being drawn exactly once --the max on the plot is around p(drawn exactly once)=.5. Jun 22, 2013 at 1:49

EDIT after seeing your code for the simulation

I think you did not consider the fact that a certain number can also occurr more than once within one sample. After accounting for that, we get the same results. Here is your modified code:

fx02<-function(ll,n,k){
a1<-matrix(0,n,2)
samp1 <- sample(1:n, k, replace=TRUE)
samp2 <- sample(1:n, k, replace=TRUE)

a1[sort(unique(samp1)), 1] <- as.numeric(table(samp1))
a1[sort(unique(samp2)), 2] <- as.numeric(table(samp2))

sum(rowSums(a1)==1)/n
}

ss<-(1:60)*5   #The grid of values of k for which we'll compute the probability.
a4<-matrix(NA,length(ss),2)
for(i in 1:length(ss)){
a2<-ss[i]
a3<-c(lapply(1:1000,fx02,n=100,k=a2),recursive=TRUE);
a4[i,]<-c(a2,mean(a3))
}
plot(a4,xlab="k",ylab="frequency of distinct draw", pch=16, las=1)
abline(h=0.3697296, v=50)

max(a4[,2])
 0.36971 For the experiment 2 (with replacement), I think the probability that a certain number is drawn exactly once is: $$P_{\text{once}}=2k(n-1)^{(2k-1)}(1/n)^{2k}$$ or $$P_{\text{once}}=2k(n-1)^{(2k-1)}n^{-2k}$$

This can be checked by a simple example where $n=3$ and $k=2$. The probability to draw a certain number exactly once is in this case $2(4/9)^{2}\approx0.395$.

The maximum occurs of the formula above occurs for: $$k_{\text{max}}=\lfloor-\frac{1}{[2(\log(n + 1) + \log(1/n))]}\rceil$$ or $$k_{\text{max}}=\lfloor-\frac{1}{[2(\log(n + 1) - \log(n))]}\rceil$$ Where $k_{\text{max}}$ is rounded to the nearest integer. This is roughly $n/2$ for larger $n$, as already mentioned in the comments by @user603. For $n=100$, $k_{\text{max}}\approx 49.75$ so 50. The maximum probability for $n=100$ and $k_{\text{max}}=50$ would then be around $0.3697$ (as already worked out in the comments by @user603).

I set up a simulation to check this result in R:

prob.once <- vector()

draw.once <- function(n, k, sim=10000, repl=TRUE){

for ( i in 1:sim ) {

samp1 <- sample(1:n, size=k, replace=repl)
samp2 <- sample(1:n, size=k, replace=repl)

if ((is.element(1, samp1) & !is.element(1, samp2) & !is.element(1, samp1[duplicated(samp1)])) |
(!is.element(1, samp1) & is.element(1, samp2) & !is.element(1, samp2[duplicated(samp2)])) ){
prob.once[i] <- 1

} else {
prob.once[i] <- 0
}
}
mean(prob.once)
}

krepl <- 1:300
probs.repl <- sapply(krepl, FUN=draw.once, n=100, sim=20000, repl=TRUE)

plot(probs.repl~krepl, pch=16, type="p", lwd=2, las=1, ylab="Probability", xlab="k", col="steelblue")

abline(h=0.3697296)
abline(v=50) The simulated result seems to confirm the above considerations.

• I can't run your code as I get an error message: "Error in prob.once[i] <- 0 : object 'prob.once' not found". Jun 22, 2013 at 12:57
• @user603: Sorry, I forgot to add the initializer for that vector. Please try again. Tell me if it doesn't work again. Jun 22, 2013 at 15:03
• @JoshO'Brien I took the first derivative of the $P_{\text{once}}$ and set it to zero and solved it. I used my calculator to do this though. The first derivative is: $2(n-1)^{2k-1}(1/n)^{2k}(1+2k(\log(n+1)+\log(1/n)))$. Sorry that I can't give you more details. Jun 22, 2013 at 20:05
• @JoshO'Brien You could use WolframAlpha to do such calculations online. It's free. Look at that. Jun 22, 2013 at 20:22
• @COOLSerdash Terrific! I had been able to work out the first derivative using the chain rule, but was then completely stuck. I guess WolframAlpha is the calculator I was wishing for! Jun 22, 2013 at 20:34