(This article was originally published at R snippets, and syndicated at StatsBlogs.)
And the question was whether to calculate it exactly or simulate the result. With GNU R you can do both. However - not surprisingly - simulation is much easier and faster to perform and exact calculation is coding error prone. Have a look at the codes.
First start with simulated result. The code is clean and simple:
set.seed(1)
deck <- rep(c(1:4, rep(0, 9)), 4)
points <- replicate(1000000, sum(sample(deck, 26)))
print(2* mean(points > 32))and we get the following result:
[1] 0.00687
On the other hand exact result requires careful counting of all possible card combinations:
library(caTools)
result <- data.frame(HCP=0:40, exact = 0)
result[1, 2] <- choose(16, 0) * choose(36, 26)
for (i in 1:16) {
HCP.sums <- table(rowSums(combs(rep(1:4,4), i))) *
choose(36, 26 - i)
levels <- as.integer(names(HCP.sums)) + 1
for (i in 1:length(HCP.sums))
result[levels[i], 2] <- result[levels[i], 2] +
HCP.sums[i]
}
result[, 2] <- result[, 2] / sum(result[, 2])
print(2* sum(result$exact[result$HCP > 32]))It took: (a) some thinking before coding, (b) 10x more time to code and (c) one mistake in the process to get the desired results. The output of the procedure is the following:
[1] 0.0069593
So the conclusion is that some pair gets a hand good enough for 6NT on the average once in 150 games.
I wanted check sure whether the simulation and exact results are similar so I decided to plot the resulting HCP distributions using the code:
par(mar=c(4, 4, 1, 1))
result$sim <- sapply(0:40, function(x) { mean(points == x) })
plot(result$HCP, result$exact, pch = 4,
xlab = "HCP", ylab = "Probability")
points(result$HCP, result$sim, pch = 3, col ="red")Tthe result given below shows that simulation approximates exact result pretty well.
Please comment on the article here: R snippets
