7
$\begingroup$

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$):

enter image description here

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")
$\endgroup$
  • 1
    $\begingroup$ Are you trying to maximize the chance of an individual ball being drawn exactly once or the number of balls being drawn exactly once? $\endgroup$ – soakley Jun 22 '13 at 0:34
  • $\begingroup$ Important point. The chance. How can I re-phrase my question to make this clear? $\endgroup$ – user603 Jun 22 '13 at 0:36
  • 3
    $\begingroup$ 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? $\endgroup$ – soakley Jun 22 '13 at 1:28
  • $\begingroup$ @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$. $\endgroup$ – Josh O'Brien Jun 22 '13 at 1:36
  • 1
    $\begingroup$ @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. $\endgroup$ – user603 Jun 22 '13 at 1:49
4
$\begingroup$

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])
[1] 0.36971

user603 code plot


Original answer

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)

Simulation plot

The simulated result seems to confirm the above considerations.

$\endgroup$
  • $\begingroup$ I can't run your code as I get an error message: "Error in prob.once[i] <- 0 : object 'prob.once' not found". $\endgroup$ – user603 Jun 22 '13 at 12:57
  • $\begingroup$ @user603: Sorry, I forgot to add the initializer for that vector. Please try again. Tell me if it doesn't work again. $\endgroup$ – COOLSerdash Jun 22 '13 at 15:03
  • 2
    $\begingroup$ @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. $\endgroup$ – COOLSerdash Jun 22 '13 at 20:05
  • 1
    $\begingroup$ @JoshO'Brien You could use WolframAlpha to do such calculations online. It's free. Look at that. $\endgroup$ – COOLSerdash Jun 22 '13 at 20:22
  • 1
    $\begingroup$ @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! $\endgroup$ – Josh O'Brien Jun 22 '13 at 20:34

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.