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

print**<-**replicate**(**1000000, sum**(**sample**(**deck, 26**)))****(**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

print**[**, 2**]****<-**result**[**, 2**]****/**sum**(**result**[**, 2**])****(**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

points**=**"HCP", ylab**=**"Probability"**)****(**result

**$**HCP, result

**$**sim, pch

**=**3, col

**=**"red"

**)**

Tthe result given below shows that simulation approximates exact result pretty well.

Not bad -- now write a version for suited slams, where you must (at the very least) throw in point counts for short (<2 cards) and long( >6 cards) suits. :-)

ReplyDeleteNice. Another extension could be to check for control of all four suits. For example, either all four Aces or three Aces plus a protected King for a small slam in No Trump. In suited slams, you'd need to include singletons and voids. Much higher level of programming difficulty, I guess.

ReplyDeleteAnother thought: you've calculated the chances of getting *at least* 6NT. Properly speaking, you should calculate the chance of getting *exactly* 6NT by subtracting the chance of getting 7NT.

ReplyDeleteI wanted a hand "good enough" for a small slam.

Deleteprint(2 * mean(points > 32))

ReplyDeleteWhy multiply by 2? Is it because you want the probability of either partnerships having as slam hand as opposed to the probability of just you and your partner have a slam hand?

Yes - I count the probability that either partnership has at least 32 points.

DeleteYou could alternatively write:

print(mean(pmax(points, 40 - points) > 32))

This is equivalent because the distribution is symmetrical and the events are mutually exclusive.