[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Recently I played bridge with my friends. Being frustrated with several consecutive poor hand distributions we asked ourselves a question what is the probability of having a hand good enough for a small slam. A well known rule of thumb is that you need 33+ HCP for 6NT. But we could not find information about the probability of such an event. So we decided to calculate it.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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)< o:p>
deck <- rep(c(1:4, rep(0, 9)), 4)< o:p>
points <- replicate(1000000, sum(sample(deck, 26)))< o:p>
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)< o:p>
result <- data.frame(HCP=0:40, exact = 0)< o:p>
result[1, 2] <- choose(16, 0) * choose(36, 26)< o:p>
for (i in 1:16) {< o:p>
HCP.sums <- table(rowSums(combs(rep(1:4,4), i))) *< o:p>
choose(36, 26 – i)< o:p>
levels <- as.integer(names(HCP.sums)) + 1< o:p>
for (i in 1:length(HCP.sums))< o:p>
result[levels[i], 2] <- result[levels[i], 2] +< o:p>
HCP.sums[i]< o:p>
}< o:p>
result[, 2] <- result[, 2] / sum(result[, 2])< o:p>
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))< o:p>
result$sim <- sapply(0:40, function(x) { mean(points == x) })< o:p>
plot(result$HCP, result$exact, pch = 4,< o:p>
xlab = “HCP”, ylab = “Probability”)< o:p>
points(result$HCP, result$sim, pch = 3, col =“red”)Tthe result given below shows that simulation approximates exact result pretty well.
To leave a comment for the author, please follow the link and comment on their blog: R snippets.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.