HCP odds April 17, 2007
Posted by dorigo in games, mathematics, personal.trackback
Many thanks to Riqie Arneberg, who criticized me for blindly trusting some information on a web site about the odds of picking up 13 cards in a 52-card deck and ending up with 28 or more high-card points (HCP). She was right: if you do not know something and your life is not at stake, go search the web; but if you can find out the solution of your problem by yourself, trust nobody…
For those of you that do not know what is a HCP count: Aces count 4 points, Kings 3, Queens 2 points, and Jacks 1. All other cards do not count. The hand I was dealt yesterday was a whooping 28 HCP, since I had 3 aces, 4 kings, and 2 queens.
So it took me ten minutes of fun this morning to write the following piece of code:
program prob_hcp
c
c Computes odds of getting a given HCP in a 13-card hand
c ——————————————————
real rndm, dummy
integer i,j
integer hcp
integer x
real sum_hcp
real total
real seen(40)
integer dealt(13)
c
c Start of code – some initializations
c ————————————
total = 0
do hcp = 1, 40
seen(hcp) = 0
enddo
c
c Loop here on hands dealt
c ————————
10 sum_hcp = 0
c
c Refresh deck
c ————
do j=1, 13
dealt(j)=0
enddo
c
c Deal hand
c ———
do i = 1, 13
c
c Deal a new card
c —————
20 x = 1+ifix(rndm(dummy)*52)
c
c Check the card is still in the deck
c ———————————–
if ( i.gt.1 ) then
do j = 1, i-1
if ( x.eq.dealt(j) ) goto 20
enddo
endif
c
c Update deck
c ———–
dealt(i)=x
c
c Add HCP of cards to sum
c ———————–
if ( x.le.4 ) then ! an Ace
sum_hcp=sum_hcp+4
else if ( x.le.8 ) then ! a King
sum_hcp=sum_hcp+3
else if ( x.le.12 ) then ! a Queen
sum_hcp=sum_hcp+2
else if ( x.le.16 ) then ! a Jack
sum_hcp=sum_hcp+1
endif
enddo
c
c Update total hands dealt
c ————————
total=total+1
c
c Print out results
c —————–
if ( mod(total,100000).eq.0 ) then
print *, ‘****** Total hands dealt = ‘, total, ‘ ******’
endif
c
c Print probability that hand has at least hcp points
c —————————————————
do hcp = 1, 40
if ( sum_hcp.ge.hcp ) seen(hcp)=seen(hcp)+1
if (mod(total,100000).eq.0 ) print *, hcp, seen(hcp)/total
enddo
goto 10
c
c That’s all folks!
c —————–
end
Here is the output after 16.7 million hands dealt:
——— k = 16700000.——–
1 0.99636215
2 0.988477051
3 0.974916101
4 0.950280249
5 0.911806643
6 0.86003077
7 0.79450506
8 0.714180052
9 0.625258327
10 0.531646669
11 0.437651783
12 0.348090291
13 0.267864734
14 0.198694006
15 0.141775325
16 0.0975567698
17 0.064420417
18 0.0407968275
19 0.0247638319
20 0.0144373057
21 0.00801041909
22 0.00424473034
23 0.00214544917
24 0.00102688628
25 0.000462155702
26 0.000193892221
27 7.59880204E-05
28 2.86826344E-05
29 9.70059864E-06
30 2.75449111E-06
31 5.38922166E-07
32 5.9880243E-08
It transpires that P(>=28HCP) = 28.68E-6, or a bit less than three times in a hundred thousand.
[...] I did the computation myself ( see the next post). I find a probability of 2.9 times in a hundred [...]
Hi Tommaso.
The site that you pointed gives the probability for each number of points while you give the probability of 28 or higher. I guess both gives the same answer as you can see by summing the last entries in their site.
Hi Pedro,
you are right – I looked at it a bit too quickly yesterday. I am actually happy to see my program is bug free…
Cheers,
T.
A comment I want to make myself is, is there an easy way to compute (even as an approximation) the HCP probability with probability calculus, without reverting, that is, to a random number generator ?
I find it hard to do it. Take P(28): you can do 28 points in the following ways:
1) 4A,4K;
4A,4Q,4J;
2) 4A,3K,1Q,1J;
3) 4A,3K,3J;
4) 4A,2K,3Q;
5) 4A,2K,2Q,2J;
6) 4A,1K,4Q,1J;
7) 4A,1K,3Q,4J;
9) 3A,4K,2Q;
10)3A,4K,1Q,2J;
11)3A,4K,4J;
12)3A,3K,3Q,1J;
13)3A,3K,2Q,3J;
14)2A,4K,4Q;
15)2A,4K,3Q,2J;
16)2A,4K,2Q,4J;
17)2A,3K,4Q,3J.
Now, we have to sum the probability of each of those 17 classes. And that, too, is not trivial…. I think my interest for this problem is rapidly dissolving
T.
Dorigo:
I would speculate that both the web page you cited, and your own calculation, are correct. There is no reason to suspect either. The problem is you mis-interpret the questions. I see the two sets of numbers are consistent with each other, if you look again exactly how the question is asked. Some one as smart as you should be able to figure it out without any further hint.
Didn’t read the previous comment. Looks like Petro already pointed it out.
Yes, I wrote about it above.
Cheers,
T.
To answer dorigo’s question, a simple recursive routine can measure the precise probability. As usual, it is simpler to answer the (aparently) more general question of which is the prob. of drawing at least p points in a hand of n cards from a deck with a aces, k kings, q queens, j jacks and x non-point cards (Here’s a VBA example (it’s not optimized at all so it takes 15min or so to run on my pc and gives an answer of 28.2644358197…E-6):
Function PR(a, k, q, j, x, p, n)
If (a 4 * n Then PR = 0
Else
t = a + k + q + j + x
PR = PR(a – 1, k, q, j, x, p – 4, n – 1) * a / t + _
PR(a, k – 1, q, j, x, p – 3, n – 1) * k / t + _
PR(a, k, q – 1, j, x, p – 2, n – 1) * q / t + _
PR(a, k, q, j – 1, x, p – 1, n – 1) * j / t + _
PR(a, k, q, j, x – 1, p, n – 1) * x / t
End If
End Function
To try for the problem at hand, simply run:
Sub test()
MsgBox (PR(4, 4, 4, 4, 36, 28, 13))
End Sub
DanC/
Sorry, the code didn’t paste right in the previous. Here it is:
Function PR(a, k, q, j, x, p, n)
If (a 4 * n Then PR = 0
Else
t = a + k + q + j + x
PR = PR(a – 1, k, q, j, x, p – 4, n – 1) * a / t + _
PR(a, k – 1, q, j, x, p – 3, n – 1) * k / t + _
PR(a, k, q – 1, j, x, p – 2, n – 1) * q / t + _
PR(a, k, q, j – 1, x, p – 1, n – 1) * j / t + _
PR(a, k, q, j, x – 1, p, n – 1) * x / t
End If
End Function
DanC/