```
rm(list=ls())
# number of simulations
n.sims <- c(10, 50, 100, 250, 500, 1000, 1500, 2000, 3000, 5000, 7000, 10000)
Est.Prob <- numeric(length(n.sims))
kids <- c(2, 3, 5, 7,
11, 13, 17, 19,
23, 29, 31, 37)
for (i in 1:length(n.sims)){
n.sim <- n.sims[i]
n.times <- numeric(n.sim)
for (j in 1:n.sim){
table.1 <- matrix(0, 10, 4)
table.2 <- matrix(0, 10, 4)
table.3 <- matrix(0, 10, 4)
for (k in 1:10){
positions <- sample(kids, replace=FALSE)
table.1[k,] <-positions[1:4]
table.2[k,] <-positions[5:8]
table.3[k,] <-positions[9:12]
} #next k
Table <- data.frame(Table.1=apply(table.1, 1, prod),
Table.2=apply(table.2, 1, prod),
Table.3 =apply(table.3, 1, prod))
n.times[j] <- sum(diff(sort(rowSums(Table)))==0)
} # next j
Est.Prob[i] <- sum(n.times)/n.sims[i]
} # next i
plot(Est.Prob~n.sims, type='o', pch=21, col='blue', bg='yellow',
ylab='Estimated Probability', xlab='Sample Size', ylim=c(0,0.05),
main='Add Title Here')
```