First read the code and try to figure out what each step is doing. Then execute the code and see results. Remember that you can select, copy, and paste R commands directly into the R session window.
k <- 4 # no. of friends
n <- 15 # no. of restaurants
x <- 1:n
m <- 10000 # no. of iterations of loop
#
y <- numeric(m)
for (i in 1:m) {
y[i] <- length(unique(sample(n,k,replace=TRUE)))
}
sum(y < k) / m
The answer is close to .353: Now, let's use R to compute exact and approximate answers for a small set of values of k:
k.values <- c(10,20:25,30,40,50,100) #table these values of k
p.n.k <- numeric(length(k.values))
for (i in 1:length(k.values)) {
v <-seq(1-1/n,1-(k.values[i]-1)/n,-1/n) #set up factors
p.n.k[i] <- 1-prod(v) #compute the p_n,k probabilities
}
approx <- 1-exp(-k.values*(k.values-1)/(2*n))
Here is a table of results, with values of k, exact answers in column
2, and approximate answers in column 3.
k exact approx 10 0.1169482 0.1159907 20 0.4114384 0.4058051 21 0.4436883 0.4374878 22 0.4756953 0.4689381 23 0.5072972 0.5000018 24 0.5383443 0.5305363 25 0.5686997 0.5604122 30 0.7063162 0.6963200 40 0.8912318 0.8819900 50 0.9703736 0.9651313 100 0.9999997 0.9999987