Playing with the pack of cards on my desk, something just happened that I didn't expect!
Look at the deck with the cards face-up. Here's a rule: look at the top card. Draw that many cards from the deck. Repeat until the deck is empty.
I did this, and dealt out exactly the right number of cards (at the end, I had 3 cards left and the top card was a 3).

What's the probability of this happening?

· · Web · 6 · 4 · 8

@christianp How much are J, Q and K worth? 11, 12 and 13? (Anyway, I thought about it a little bit and couldn’t find how to control the combinatory explosion.)

@christianp That's single player blackjack with 52 instead of 21, innit?

@christianp Since the deck is shuffled you can transpose the "look at the next card" and "remove the required number of cards" steps. This means it's equivalent to drawing cards from the top until they either sum to 52 (success) or higher (failure).

This doesn't really give a full answer, only simplifies the problem a little.

@christianp Could it be a lot of annoying logic to end up at 1:13?

@LovesTha @christianp My guess would be 1/e, but that's just what I always guess when it's a problem about finding one exact match when a lot of possible matches might but don't have to exist.

@christianp If each card appears randomly with replacement (magically?) and you reach a pile with k cards left, let P_k = probability you won. P_k = 0 if k < 0 (you ran out of cards in last draw), P_0 = 1 (you won!), and otherwise P_k = sum_{i=1}^13 P_{k-i}/13. A quick computational experiment shows that for large k like k=52, P_k appears to converge to 1/7.

@11011110 @christianp My simulation came up a little bit short of 1/7, so this looks good to me.

@11011110 @christianp Following @timorl's logic, I think we only need to check (at most) a few trillion cases -- I think once you have dealt 18 cards, you've either reached 52 or know you can't hit it.

@icecolbeveridge @11011110 @christianp I did a quick numerical simulation this morning too and came up with .1422, which made me think, "Wait, isn't that 1/e^2?" It's not, but you can see 1/e^2 from there.

@nebusj @11011110 @christianp

OK, a bit of thought that I've not yet coded up: there are <10^6 integer partitions of 52, and only ~42,000 with no more than four matches and a max <= 13.

Loop through those, summing (perms of partition)\(\times\)(perms of remainder). Divide the result by 52! and I think we're done.

@nebusj @11011110 @christianp Oops, off by a factor of \( 24^{13} \), no biggie.

I get 0.1420342593977892,

@icecolbeveridge @nebusj @11011110 nice! So if I did this with a class of kids, I'd expect to see it happen a few times

@christianp @icecolbeveridge @nebusj Here's a conceptual argument why the limit (for drawing with replacement) is exactly 1/7. If you travel along the number line, drawing a card to determine each step size, then the average step length is \(\sum_{i=1}^{13}\tfrac{i}{13}=7\). By Chernoff, whp traveling \(N\) steps covers distance\(7N+o(N)\), and hits an average of 1/7 of the numbers along the way. By the time you near an empty deck, randomization will have ironed out fluctuations in this average.

@11011110 @christianp @nebusj I like it. Rephrasing in terms I can understand: the cards range from 7-6 to 7+6, so on average each card moves you a 7-count forward and the variation from that spreads out evenly. (I acknowledge the hand-waving in that.)

@icecolbeveridge @christianp @nebusj Chernoff bounds are just a technical way of saying that, if you make a lot of small random steps, the average of the sizes of the steps you actually took is very likely to be close to its expected value.

The "ironed out fluctuations" part can be made rigorous (but a lot more technical) by describing it as converging to the stable distribution of a Markov chain on a de Bruijn graph, whose states describe which numbers get hit within a 13-number window.

@11011110 @christianp

Quick question: how do you consider the effect that you only have 4 suits for every card?

@sojournTime @christianp I'm not considering it. That's what I meant about "with replacement". It's an assumption that makes the mathematics easier but also less accurate for the actual problem.


So I've been thinking about the problem almost like a knapsack problem. Let \(\vec{c}=[1,2,3,4,5,6,7,8,9,10,11,12,13]^T\), which is the vector of card values, and \(\vec{a}=[a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8,a_9,a_{10},a_{11},a_{12},a_{13}]^T\) where \(\forall_{i \in \{1,\dots,13\}} a_i \in \{0,\dots,4\}\), which are the number of individual cards (e.g., 4 aces, 1 twos, etc.). Then, solutions to the puzzle could be written as the following:

\(\min (\vec{c}^T\vec{a}-52)^2\)


\(\vec{a}^T\textbf{1} \leq 52\)

Thus, any solution found by this minimization problem would have \(\vec{c}^T\vec{a}=52\). Then, given a solution, the number of combinations for the solution is as follows:

\((\vec{a}^T\textbf{1})! \times (52-\vec{a}^T\textbf{1})! \times \prod_{i \in \{1,\dots,13\}}(4 nCr a_i)\)

The total number of combinations is the sum of all of these values.

Of course, solving the minimization problem is difficult, but I would assume just like the linear knapsack problem, dynamic programming, as @narain mentioned before, would do the trick. Please let me know if I have made any mistakes.

@christianp @narain

Alright, so I couldn't help myself and ran the optimization problem in AMPL 😅 .

Before I get to my results, I want to make one more observation I made. Since we are dividing all the combinations by 52!, the solutions follow a multivariate hypergeometric distribution:

\(\frac{\prod_{i \in \{1,\dots,13\}}(\text{4 nCr} a_i)}{(\text{52 nCr }\vec{a}^T\textbf{1})}\)

Now, when I ran the optimization algorithm through AMPL, I got a total of 41,846 solutions of the form \(\vec{a}=[a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8,a_9,a_{10},a_{11},a_{12},a_{13}]^T\). I ported over my answers into a csv file so that I can do the combination/probability calculations in Python. Just for my own verification, I checked to see if all of the solutions had a sum of 52 and were all unique combinations (my trust in getting the right settings in AMPL is not high 😅 ), which thankfully they were. Then, I calculated the probabilities by the number of combinations divided by 52! and by the multivariate hypergeometric distribution method. Then, the sum of each of these probabilities resulted in the total probability.

So, the answer I got was...


This number seems very high compared to everyone else's simulations, so if anyone wants to see my code, please let me know, and I will post it to GitHub. I have attached a screenshot of my Python results that includes the total number of combinations and the probabilities. With that said, I really enjoyed this problem and had a lot of fun solving it!

@christianp @narain

Alright, so I went over my code with a colleague, and he found an error in my code (forgot to round my floating number before converting to integer). My original number of solutions (41,846) still holds, but my probability and combination calculations were off. My current probability now matches the simulation number of approximately 14.2%. I have uploaded my code and solution on GitHub here:

The following is an updated screenshot of my solution:

Sign in to participate in the conversation

The social network of the future: No ads, no corporate surveillance, ethical design, and decentralization! Own your data with Mastodon!