I have started a new book called The Art of R Programming by Norman Matloff and I’m really digging it. I won’t blog about each chapter the way I did Machine Learning for Hackers, but I did come across something I thought made good blog material.

In Chapter 8 (Doing Math and Simulations in R), Matloff presents an example on combinatorial simulation. The motivation for the example is the following problem:

Three committees, of size 3, 4 and 5, are chosen from 20 people. What is the probability that persons A and B are chosen for the same committee?

Says Matloff: “This problem is not hard to solve analytically, but we may wish to check our solution using simulation…” I suppose it’s not that hard, but if you don’t think about these kinds of problems on a regular basis they can trip you up. He doesn’t give the analytic solution but does provide code for a neat simulation of it:

# 8.6.3 combinatorial simulation

sim <- function(nreps) {commdata <- list()commdata$countabsamecomm <- 0for (rep in 1:nreps) {commdata$whosleft <- 1:20commdata$numabchosen <- 0commdata <- choosecomm(commdata,5)if (commdata$numabchosen > 0) nextcommdata <- choosecomm(commdata,4)if (commdata$numabchosen > 0) nextcommdata <- choosecomm(commdata,3)}print(commdata$countabsamecomm/nreps)
}

choosecomm <- function(comdat,comsize) {committee <- sample(comdat$whosleft, comsize)comdat$numabchosen <- length(intersect(1:2,committee))if (comdat$numabchosen == 2)comdat$countabsamecomm <- comdat$countabsamecomm + 1comdat$whosleft <- setdiff(comdat$whosleft,committee)return(comdat)
}

That’s quite a bit of code. The first block is the main function, called “sim”. The second block is the other function, “choosecomm”, that gets called in the first function. The whole thing makes creative use of a list and the intersect() and setdiff() functions. If you want to simulate selecting three committees of 3, 4 and 5 people 10,000 times and see how many time persons 1 and 2 are on the same committee, enter sim(100000) at the command line. You should get about 0.10.

Now how do we solve this using math and probability? First off, instead of 20 people, I like to think of 20 chips, 18 of which are white and 2 that are red. Imagine they’re in a bag and we reach in and scoop out 12 and immediately and randomly divide into three groups of 3, 4 and 5. What is the probability one of those groups has both red chips? To calculate this we need to enumerate the possible selections that include both red chips and divide that by all possible selections . First let’s do it for the committee of 5. All possible combination of 5 from 20 is . All possible combinations including both red chips are . We can calculate this in R using the choose() function as follows:

c5 <- choose(18,3)/choose(20,5)

We then repeat this for the other two committees of 4 and 3, like so:

c4 <- choose(18,2)/choose(20,4)
c3 <- choose(18,1)/choose(20,3)

Summing this up we get c3 + c4 + c5 = 0.10. Now it may seem strange that we always set up our denominator as if we’re drawing from 20. After all, in the simulation above, there’s a specific order. First 5 are chosen. If the two people are not both on that committee, then we draw 4 and then 3. It seems that should be taken into account when you solve this mathematically. But you don’t have to. That’s why I like to imagine just scooping 12 people up and assigning instant membership to them on a random basis. While that’s not how committee selection usually happens in “real life” it makes thinking about an analytic solution easier.

引用: http://www.clayford.net/statistics/a-combinatorial-simulation/