The following problem is from the Fiddler (the “spiritual successor to FiveThirtyEight’s The Riddler column”), involving the game of tenpin bowling with traditional scoring:

If you knock down 100 pins over the course of a game, you are guaranteed to have a score that’s at least 100. But what’s the minimum total number of pins you need to knock down such that you can attain a score of at least 100?

I thought this was an interesting problem, in part because traditional scoring in bowling is a bit of a mess. Every pin that you knock down scores at least one point, but some pins may yield additional bonus points… based on previously knocked down pins.

My first thought was to transform the puzzle into an integer linear programming (ILP) problem. In hindsight, the solution ends up having a nice enough form to suggest a more elegant cocktail-napkin proof. But those scoring rules are just messy enough to want some more brute-force verification from the computer. Or at least that’s what I told myself to justify the interesting coding problem.

The Mathematica code is on GitHub; it’s not much, just a few dozen lines. The solution is simple to describe: we can score at least 100 points by knocking down just 44 pins, rolling four consecutive strikes, immediately followed by knocking down four pins, and no more, for a total score of 30+30+24+14+4=102.

We certainly can’t fully brute-force the problem by examining all possible games: the number of distinct “scoreboards” in a game of bowling is

However, at the other extreme, I also wasn’t able to reduce the problem to just a single ILP instance. Not only is the total score (which we are bounding) not a linear function of the variable numbers of pins knocked down by each roll, but even the total number of pins knocked down (which we are trying to minimize) is not a simple linear function of the variables, because of the possible extra frames after a strike or spare in the tenth frame.

Instead, we can trade execution time for preserved code simplicity, by evaluating not just one ILP instance, but of them, by grouping game outcomes according to the sequence of frame “types”: strike, spare, or open. For example, the eventual solution described above occurs in a game of the form (strike, strike, strike, strike, open, open, open, open, open, open). Having fixed the type of each frame in the game, now we can describe the number of pins knocked down and the resulting total score as linear functions of variables indicating the number of pins knocked down by each roll.

Introduction

Let’s play a game: I will repeatedly flip a fair coin, showing you the result of each flip, until you say to stop, at which point you win an amount equal to the fraction of observed flips that were heads. What is your strategy for deciding when to stop?

This weekend 6/28 is “Two-Pi Day” (or Tau Day, if you like). Celebration in the classroom can take many forms: for example, we can write a program to calculate an approximation of the value of using Monte Carlo simulation. Or we can find some mathematical problem where the occurrence of is surprising or unexpected. I think the game described above is an interesting case where we can do both at the same time.

Throwing darts

The most common approach to estimating using Monte Carlo simulation involves selecting random points uniformly distributed in a square; imagine throwing darts at a square dartboard, as shown in the figure below.

Let the indicator random variable equal 1 if a dart lands within the inscribed circle– or within the inscribed 90-degree circular sector shown above, which has the same area– otherwise . Then is the ratio of areas , so that we can estimate by Monte Carlo estimation of , where . In Python:

import random

def dart():
    x, y = random.random(), random.random()
    return 4 * int(x * x + y * y < 1)

n = 1000000
print(sum(dart() for _ in range(n)) / n)

Note how simple the simulation of dart() is, just a couple of random numbers and a few arithmetic operations. Indeed, this might seem perhaps too simple: we don’t really need the randomness of Monte Carlo simulation here, and could instead integrate over an evenly spaced grid of thrown darts. And despite the simple arithmetic implementation of the simulation, the “pi-ness” of the process– there is a circle clearly visible in the diagram– is rather unsurprising.

Throwing noodles

Another common approach to estimating is Buffon’s noodle. Imagine throwing uncooked straight spaghetti noodles, each one foot long, onto a floor marked with parallel lines spaced one foot apart, with the random variable equal to the number of times a particular noodle crosses a line on the floor.

It turns out that . (Note this is an expected value, not just a probability, because it works even if we cook the noodles, so that can take on non-negative integer values beyond just 0 and 1.) We can similarly estimate by Monte Carlo simulation of this noodle-throwing process, appropriately transforming the resulting estimate of . Again in Python:

import random
import math

def noodle():
    y = random.random()
    angle = math.pi * random.random()
    angle_0 = math.asin(y)
    return int(angle_0 < angle < math.pi - angle_0)

n = 1000000
print(2 / (sum(noodle() for _ in range(n)) / n))

As a method of computing an approximation of , this feels a bit circular (!) to me. We are using math.pi, the very value we are trying to estimate, in the implementation of the simulation! (I’m also skipping over the error analysis, which is more complicated than for the other simulation approaches discussed here, due to the reciprocal in the estimate.)

Flipping coins

Coming back now to the Chow-Robbins coin-flipping game, unfortunately fully specifying an optimal strategy with maximum possible expected return (known to be approximately 0.793) remains a complex open problem. But let’s consider a strategy that, although it is suboptimal, is very simple to describe, still performs reasonably well, and turns out to have interesting behavior relevant to this discussion: let’s suppose that we stop when the number of heads first exceeds the number of tails , so that we win .

Then it’s a nice exercise to show that our expected return using this strategy is

So we can estimate in the same way we did with the dartboard. In Python:

import random

def coin_flips():
    heads = 0
    flips = 0
    while heads <= flips - heads:
        heads += random.getrandbits(1)
        flips += 1
    return 4 * heads / flips

total = 0
n = 0
while True:
    total += coin_flips()
    n += 1
    print(f'Estimate after {n} iterations = {total / n}')

This works… sort of. There are a lot of interesting things happening here. First, the variance of our estimate is less than one-third that for a dart throw, so that we should need correspondingly fewer Monte Carlo iterations to achieve comparable accuracy.

Also, each individual coin flip is a lot less computationally expensive than each dart throw. The latter requires 128 random bits to generate two floating-point values for the position of a dart, while a coin flip is a single random bit. The ability and relative speed of extraction of a single random bit varies from language to language and library to library. For example, even sticking to Python, it’s interesting to compare the above random.getrandbits(1) with the theoretically equivalent but much slower random.randint(0, 1).

But these potential benefits are dwarfed by the overall cost of each Monte Carlo iteration simulating a single playout of the game, which generally involves more than just a single coin flip. You may notice that the output from the above code scrolls by with not-really-blazing speed… but periodically there is a significant pause, sometimes for a long time.

The problem is that the expected number of coin flips– that is, the average duration of a single Monte Carlo iteration– is infinite. In other words, even if we know that we need, say, roughly 12 million iterations of the game to achieve 3 digits of accuracy in our estimate of , we can’t even estimate how long those 12 million iterations will take to execute.

All is not entirely lost, however. Note that although the expected duration of a game is infinite, the program above is never in any danger of being well and truly hung in the inner loop. That is, even during one of those long pauses, we are guaranteed— with probability 1– to observe eventual completion of the current game and subsequent starting of the next new game.

We are right on the razor’s edge of hanging, though. That is, suppose that we play the game with a biased coin, so that the probability of heads is . Then there is a positive probability that a game will never end, with a meandering deficit of heads that is never overcome. (Which means that the endless loop through games above is guaranteed, with probability 1, to eventually encounter and get stuck on such a never-ending game.)

On the other hand, if , then the situation improves. For example, let’s suppose that we play the game with a biased coin with . Then the expected return from this game is , so we can still estimate by playing this modified game:

import random
import math

def coin_flips():
    heads = 0
    flips = 0
    while heads <= flips - heads:
        heads += int(random.getrandbits(2) != 0)
        flips += 1
    return 2 * math.sqrt(3) * heads / flips

n = 120000
print(sum(coin_flips() for _ in range(n)) / n)

Better yet, not only is the expected number of coin flips (i.e., random bits) per game finite, it’s just 2 flips on average per game, so that each game is relatively fast to simulate. Better still, the variance is now less than 12% of the dartboard’s, so that we don’t even have to play as many games to achieve the same accuracy.

Introduction

This is a follow-up to an article from nearly a decade ago, answering what I didn’t recognize at the time as an interesting question, buried in an otherwise offhand remark. To repeat the setup, suppose that we shuffle together decks of playing cards, and play a series of rounds of blackjack, each using the same fixed CDZ- playing strategy, stopping at the “cut card”, a marker inserted in the shoe indicating when to reshuffle. Let’s place the cut card at “penetration” , that is, stop and reshuffle when at most decks or 78 cards remain in the shoe. Let the random variable be the outcome of the n-th round. The figure below shows estimates of , in percent of initial wager, from just 1000 simulated shuffled shoes.

That steep drop in expected return late in the shoe is the “cut card effect.” An intuitive, hand-waving explanation goes something like this: although the position of the cut card is fixed, the number of rounds needed to reach it varies from one shuffled shoe to the next. If near the end of the shoe we have played a larger-than-typical number of rounds, then each of those rounds must have consumed fewer-than-typical number of cards on average, and so were likely rich in tens. The remainder is thus poor in tens, a disadvantage for the player.

Inspection of the above figure suggests that the cut card effect “begins” at approximately round or so (for this particular choice of rules, playing strategy, and penetration). The motivation for this post is to show that the cut card effect actually “begins” much earlier in the shoe, at round in this case.

Extended True Count Theorem

The cut card effect is interesting because, at first glance, it seems to contradict a very useful theorem, first proved by Thorp, that is easy to state incorrectly, so for now I’ll state it vaguely:

Theorem: For each round played with the same fixed strategy, , provided we do not run out of cards.

In other words, compute the exact expected return from playing our fixed strategy off the top of the shoe: this is . Now compute the expected return from playing a second round, after having played a first round off the top of the shoe: this is , and the above theorem states that these two expected values are equal. Similarly, the expected return from playing the third round into the shoe is , etc.

Looking again at the figure above, the small variations in estimated expected return for small are purely due to sampling error. If we could in principle evaluate not just 1000 randomly shuffled arrangements of cards in the shoe, but all possible arrangements, computing each expected return exactly, then the above theorem states that this curve would be exactly constant for those small .

But clearly something different is happening for , so what counts as “small” ?

Integer linear programming

This is where the vague “not running out of cards” comes in. There are two conditions that must be satisfied to guarantee : we must be guaranteed to reach round , and we must be guaranteed to be able to complete round .

I think the latter condition is simpler to describe. We need to ensure that, no matter how the previous rounds play out, we have enough cards left in the shoe to complete the subsequent round , with neither the player nor the dealer running out of cards. For these rules and CDZ- playing strategy, we will consume at most 33 cards in any single round, so as long as we don’t start the round before the cut card, there is no danger of “running out of cards.”

It’s the former condition that is trickier, and is the heart of the problem here: we must, with probability 1, complete prior rounds before reaching the cut card. Otherwise, to even define the random variable , we must condition the population of possible shuffled arrangements of cards in the shoe to exclude those cases where we reach the cut card too soon.

So, how many rounds are we guaranteed to be able to play before reaching the cut card? This is an integer linear programming problem: for each of the 3,054,067 possible subsets of card ranks consumed in a round, let be the number of cards of rank consumed. Then the problem is to find non-negative integers :

minimizing

subject to

The code is on GitHub: C++ to compute all possible subsets of cards comprising the player hands in a round, and Mathematica to do the same for the dealer, to compute the outer product of player and dealer outcomes, and to solve the linear programming problem.

Interpreting results

The resulting minimum value of the objective function is 18. This means that it is possible– albeit unlikely– to reach the cut card after having played just 18 rounds. If we could eliminate the sampling error in the above figure, the actual exact expected return would be constant for .

But the extended true count theorem has nothing to say about the expected return from round or later. Empirically, these mid-shoe departures from are almost certainly very small… but they are just as almost certainly non-zero.

At a blackjack table, imagine next to the dealer is a display indicating the number of the round to be played, that increments by one after each round, and resets to 1 after each reshuffle. Further imagine that we walk up to the table in the middle of the shoe, and observe the displayed round number, with no information about any cards that have been dealt since the last shuffle. If we observe , then we can be confident that the expected return from playing the upcoming round is the same as if we walked up to a full shoe: . On the other hand, if we observe, say, , then the expected return is almost certainly not exactly equal to … but we can’t say for sure, let alone the magnitude or even the sign of the difference.

In practice, such a display does not exist; if we walk up to a table cold we don’t know how many rounds have been played up to that point. But suppose instead that we observe not the number of the round, but the number of cards left in the shoe before the cut card. How can we tell when we are in the “safe zone” of no cut card effect, with expected return exactly equal to ?

If there are 234 cards left in the shoe, before the 75% penetration cut card, that’s a full shoe, so the expected return from the round is . But suppose that there are as many as 162 cards still left in the shoe. Seems like that should still be plenty… but with even that many cards remaining, we are already “in the cut card effect,” where the extended true count theorem provides no information about the expected return from the next round. (It’s an exercise for the reader to show that this minimum of 162 cards remaining is the solution to another linear programming problem, but one that we can solve in our head.)

For the past six years, I have been maintaining (on GitHub) a machine-readable record of outcomes of the NCAA men’s basketball tournaments, over the now four decades since the current 64-team format began in 1985. (I continue to refuse to acknowledge the four Tuesday-Wednesday play-in games.) This effort was originally motivated by the annual question of the probability of picking a “perfect bracket,” i.e., correctly guessing the winners of all 63 games in all six rounds of the tournament. Despite some close calls over the years, no one has ever verifiably done this, and it is unlikely that anyone ever will.

This past post describes a method of estimating this probability for any given bracket, so that, for example, a “chalk” bracket (where a higher seeded team is always selected to beat a lower seed) is a much more likely overall outcome– albeit still unlikely as a specific outcome– than a bracket picking, say, all of the #1 seeds to be beaten by #16 seeds in the first round.

We can also apply this method to historical data, aggregating and weighting each year’s games and upsets into a single number, that we can use to compare the overall prior (un)likelihood of tournament outcomes across years, as shown in the following figure.

The black line at the bottom is , the often-quoted “one in 9.2 quintillion” probability of correctly guessing all 63 games by simply flipping a coin. The blue and red lines at the top are probabilities of chalk brackets, using two different models of individual game probabilities as a function of difference in “strength” of each team. As discussed in more detail here, blue indicates a strength following a normal distribution density, and red indicates a simpler linear strength function.

From the above figure, we see that this year’s tournament was the second most likely ever in the history of the current format. This makes sense: the Final Four teams were all #1 seeds (only the second time this has happened– the other was in 2008 with champion Kansas), and there were only 11 upsets– the fewest ever, as shown in the figure below– of a higher seed losing to a lower seed. The only more likely prior probability of overall tournament outcome was in 2007, with a then-fewest 12 total upsets, and #1 seeds Ohio State and champion Florida beating #2 seeds Georgetown and UCLA, respectively, in the Final Four.

So, how (un)likely is it that someone will eventually pick a verifiable perfect bracket? The probability of a chalk bracket– the mode of the distribution– being correct is one in roughly 100 to 200 billion. Even then, there are eight different chalk brackets… and an astronomically larger number of brackets with the dozen or so upsets that we expect to see in a given year.

I recently saw a WOW board, which seems to be a popular approach to rewarding positive behaviors in elementary classrooms. The idea is pretty simple: start with an grid (usually ) of empty squares. When a student exhibits some positive behavior, academic achievement, etc., they write their name or otherwise fill in a randomly chosen square. When the grid is “filled,” something cool happens, such as a prize drawing, extra recess, gold watches for everyone, whatever.

In most descriptions of WOW boards that I found online, the students must fill the entire grid of all 100 squares. But this particular variant was more mathematically interesting: the grid squares were numbered 1 to 100, and the teacher kept a container of slips also numbered 1 to 100; each time a student filled in a square, they didn’t simply pick any open spot on the board, but instead drew one of the slips from the container (without replacement) to determine which square to write their name in. The cool-something happens when the students first complete any row or column of the grid.

How long, on average, does this take? That is, what is the expected value of , the number of slips drawn until first completing a row or column of the grid?

This is essentially a special case of the game of Bingo discussed here before. However, the analysis approaches presented there required computing sums with an exponential number of terms: for the most brute-force approach, which we knocked down to . This is manageable for Bingo cards with , but is already unpleasant for a WOW board with .

It’s a nice problem to show that we can more efficiently compute the probability distribution– and thus the expected value– of in polynomial time using the formula

The resulting cumulative distribution for a 10×10 grid is shown below in blue, with the expected number of draws– about 71.6– shown in red.

An open question that I don’t know how to answer: what is the value of

That is, as the WOW board grid size grows large, do we expect to need to fill closer and closer to 100% of the squares before first filling a row or column?

I suspect the answer is yes, but I’m not sure how to prove it.

Introduction

I think the game of blackjack would be relatively boring, mathematically speaking, if not for splitting pairs. But that one additional playing strategy option makes analyzing the game much more challenging and interesting: if you are dealt a pair of cards with the same rank, should you “split” the pair– at a cost of another wager– and play each single-card hand separately? If you are dealt yet another pair card, should you “re-split” to a third, or even fourth, additional hand?

When evaluating splitting strategy, we usually trade some optimality for speed and simplicity. That is, we can compute the exact expected return from splitting a pair very quickly— in milliseconds– if we are willing to constrain the player’s strategy to depend only on the cards in the current hand (and maybe some limited additional contextual information), not all of the cards observed in the entire round.

But there are some analysis questions where it would be useful to know the speed of light: how much can you possibly win with truly perfect play, with every strategy decision accounting for every card dealt from the shoe? The motivation for this post is to capture my notes on an algorithm– code is on GitHub— to compute this optimal pair-splitting strategy and corresponding maximum expected return. See [2] and [3] for related work; the objective here is to significantly reduce execution time and extend the space of feasibly computable resplits.

More retrograde analysis

I was able to borrow a lot of code from the recent solution of the children’s game Gobblet Gobblers. The idea is the same here: we compute the expected value of a pair split in two steps:

1. Breadth-first search the decision tree of all possible states of the round.
2. Traverse the tree backward from leaves to the root, in order of decreasing depth, computing the expected value of each state in terms of its children.

In the first step, the breadth-first search trades the memory required to store all possible states for the speed of not repeatedly re-visiting states as in a recursive (effectively depth-first) traversal. Storing all of those states requires some careful packing; more on this shortly.

In the second step, we compute the expected value of each game state in terms of its children. For example, the value of hitting (or doubling down) from a given state is a linear combination of the values of the states reached from each possible drawn card, weighted by their probabilities; the value of a decision state is the maximum of the values obtained by standing vs. hitting, etc. By “solving” for each state in order of decreasing depth, we ensure that we have already computed the expected values of any dependent child states. We start at the bottom, with leaf states, corresponding to the player having completed all split hands in the round, which we “solve” in terms of the probabilities of outcomes of the dealer’s hand, as discussed here before.

Storing packed game states

The most interesting part of this problem was the memory management. As with Gobblet, we need a hash map from each game state to its corresponding (expected) value. And once again, we’re using the lean and mean “mask, step, index” MSI hash map implementation described by Chris Wellons [1]. (Thanks to Chris for very helpful discussions about other optimizations as well!)

But now we need even more states– well over half a billion depending on the rules and number of decks– and even more memory per state. Following is the 16-byte representation of a single game state:

struct State
{
    std::int8_t cards[11];
    std::int8_t hands[4];
    std::int8_t current;

    bool operator==(const State& rhs) const
    {
        return std::memcmp(this, &rhs, sizeof(State)) == 0;
    }
};

The cards[1] through cards[10] indicate the number of cards of each rank observed during the round, including the dealer’s up card and the initial two pair cards being split. The hands[current] indicates the current hand total, negated to indicate a soft total. The cards[0] indicates the number of cards in the current hand: either 1 for a single split pair card, 2 for a two-card hand, or 3 for more than two cards.

The hands[0] through hands[current-1] indicate totals for hands that are already completed (via stand, double down, or bust), but with a slightly different convention than hands[current]: totals are “clamped” to 16 for any total less than 17, or 17, 18, 19, 20, or 21; or 22 for any busted hand, with negation indicating that the hand was doubled down. Furthermore, this subarray of completed hand values is sorted.

The result of all of this clamping and sorting is that we cannot in general reconstruct from a State the chronological “history” of decisions and composition of cards drawn to each hand in a round. But that’s okay: the objective is to preserve only information that is exploitable by strategy decisions. When trying to determine what to do next, we don’t care whether that already-completed hand total is hard or soft, but we do care whether we doubled our wager on it, etc.

But we’re not done compressing just yet. Each element of the MSI hash map array contains not just a game state “key,” but also its expected value:

struct Solved
{
    State state;
    union
    {
        Solved* next;
        double value;
    };
};

I wrote my first computer program over four decades ago, and this is the first time I’ve found myself using a union in C++. During the first breadth-first search exploration of all possible game states, we use these 8 bytes as linked list pointers to group states by depth. Then in the second step we traverse those linked lists of states, overwriting the list pointers as we compute each state’s expected value.

One interesting thing that I learned from this project is that just zeroing memory takes a long time (“only” gigabytes per second). For this reason, you can input the desired size of the hash map, to speed up scripted evaluation of large numbers of “fast” scenarios that don’t require a large number of states. For example, to evaluate splitting 2-2 against a dealer 4 in 1D, H17, DOA, DAS, SPL1, we only need a couple of million states:

Memory allocation exponent (29=12GB, 30=24GB, etc.): 21
European no hole card rule (1=yes, 0=no): 0
Dealer hits soft 17 (1=yes, 0=no): 1
Double down after split (1=yes, 0=no): 1
Double down on any two cards (1=yes, 0=no): 1
Maximum number of splits (3=SPL3): 1
Shoe (ace through ten, -1 to exit): 4 4 4 4 4 4 4 4 4 16
Pair card: 2
Dealer up card: 4

Split 2 vs. 4:
    Searched 1247392 states to depth 17, 59% capacity.
    Solving depth 0, 100% complete.
    Elapsed time 2.675 seconds.
E(split)= 0.12079971803625923

References:

1. Wellons, C., The quick and practical “MSI” hash table [nullprogram.com]
2. ChemMeister iCountNTrack, Composition Dependent Combinatorial Analyzer [moderator on blackjacktheforum]
3. kc, Optimal Expected Values for a single split [bjstrat.net]

Suppose that couples want to participate in a Secret Santa gift exchange, where each person is secretly assigned to give a gift to one other person in the group. We need a procedure to randomize the gift-giving assignments, subject to some desired constraints:

1. Every person receives one gift; that is, the function mapping each person to their assigned gift recipient is a bijection.
2. No person is assigned to give a gift to themselves; that is, the function is a derangement.
3. No person is assigned to give a gift to their partner.
4. Every assignment generated by the procedure that satisfies the above constraints is equally likely.

A common approach to generating a random assignment is for each person to write their name on a slip of paper, put all of the slips into a hat, and pass the hat around, with each person drawing a slip to determine their assigned gift recipient. We have discussed this approach before: for example, if a person “peeks” at their slip as it is drawn and sees their own name– violating constraint (2)– then it might be tempting to “fix” the procedure inline by returning the slip to the hat and drawing another. But this turns out to violate constraint (4).

We can preserve constraint (4) by rejection sampling: if someone draws their own name– or that of their partner– then everyone immediately returns their slips to the hat, and the drawing starts again from the beginning. We can compute the probability of “accepting” such a drawing, with the added wrinkle of constraint (3) addressed here.

This post is motivated by a slightly different procedure: suppose that we have a moderator who uniformly randomly permutes the list of participants, and each person is assigned to give a gift to the person immediately following them in the shuffled list (viewed cyclically, so that the last person in the list gives a gift to the first person in the list).

This procedure satisfies constraints (1), (2), and (4), without any need for rejection re-sampling. However, what is the probability that this procedure also satisfies constraint (3)?

It’s a nice problem to compute this probability exactly; but it’s also interesting to show that this probability approaches as grows large.

Introduction

Let’s play a card game for money. For example, I (the “house”) will thoroughly shuffle together six decks of playing cards, and deal a round of blackjack against you (the “player”). That is the scenario motivating this post, but the details of the rules of the game don’t matter much for our purpose here, so we can alternatively consider a much simpler game: I will shuffle a single deck of cards, and deal one card to you and one card to me. High card wins (aces are low, face cards have value 10), with a tie going to the house.

Because we are playing for money, both the player and the house have an interest in knowing the exact (dis)advantage in the game. We can write a program to compute the player’s expected return, as a function of the composition of ranks of face-down cards in the shuffled “shoe.” For example, we might write the following Python code to compute the expected return for the high card game:

from fractions import Fraction
import numpy

def expected_value(shoe):
    """Expected value of high card game as fraction of initial wager."""
    ev = 0
    for card1 in range(10):
        for card2 in range(10):
            s = list(shoe)
            p = Fraction(s[card1]) / sum(s)
            p = p * s[card2] / sum(s)
            ev = ev + p * (1 if card1 > card2 else -1)
    return ev

print(expected_value((4, 4, 4, 4, 4, 4, 4, 4, 4, 16)))

The result is -25/169, or about 14.8% of the player’s wager going to the house.

Or is it? There is an error in the above code. Let’s leave it as an exercise for the reader to find the bug (hint: one line is missing) … and instead prove that there must be a bug, even without identifying nor fixing it.

The Extended True Count Theorem

The idea is to suppose that we shuffle the deck as usual, but then “burn” the top card of the deck, setting it aside face down and unused, before starting the game by dealing the player’s and dealer’s “hands.” Intuitively, this action should not affect the player’s expected return: neither the player nor the dealer sees the burn card, and at least in this simple high card game it wouldn’t matter even if they did, since there is no “playing strategy” that might be influenced by knowledge of the burn card’s value.

We can verify this by computing the average expected return played from the resulting 51-card deck, weighted by the probability of each possible burned top card:

def average_ev(shoe):
    """Average expected value after 'burning' the top card."""
    ev = 0
    for card in range(10):
        s = list(shoe)
        p = Fraction(s[card]) / sum(s)
        s[card] = s[card] - 1
        ev = ev + p * expected_value(s)
    return ev

print(average_ev((4, 4, 4, 4, 4, 4, 4, 4, 4, 16)))

The resulting average expected return is -557/3757: still a roughly 14.8% house edge, but not exactly the same -25/169 that we should expect, thus proving that our original implementation of expected_value(shoe) must be in error.

This approach to verifying the correctness of expected_value(shoe) applies much more generally, as described in this paper, as well as here and here. For example, instead of just burning the top card, modify average_ev(shoe) to play the game twice, computing the average expected return from the second round (effectively burning the top two cards instead of just one). Or deal cards from the top of the deck, stopping when you have either dealt an ace or a maximum of ten cards, then play the game from the resulting depleted deck. Or repeatedly flip a coin, dealing a card for each heads, stopping when you flip tails or have flipped heads a maximum of ten times, then play the game. In each case, if we ever find that average_ev(shoe) != expected_value(shoe), then we know that something is wrong.

(The converse, however, is not true: matching expected values do not imply that the implementation of expected_value(shoe) is necessarily correct. For example, let’s change the game slightly so that ties are a push instead of going to the house, with the following change to expected_value (but otherwise leaving the original bug as-is):

            ev = ev + p * numpy.sign(card1 - card2)

The reader can verify that the discrepancy goes away, with average_ev(shoe) == expected_value(shoe) == 0: not only do the values match, but they are correct, despite the bug remaining in the implementation.)

Application to blackjack

Coming back now to the original motivation for this post, let’s consider a round of blackjack with some specified rules (6D, S17, DOA, DAS, SPL3, LS), played using CDZ- strategy optimized for the full shoe. Using my blackjack software (on GitHub), we can compute the corresponding expected return from a round played off the top of the shoe (think expected_value(shoe)), as well as the average expected return after removing a single card, but playing the same fixed “full-shoe” strategy (think average_ev(shoe)):

#include "blackjack.h"
#include <iostream>

int main() {
    BJShoe shoe{6};
    BJRules rules;
    BJStrategy cdz;
    BJProgress progress;
    BJPlayer *basic = new BJPlayer{shoe, rules, cdz, progress};
    std::cout << basic->getValue() << std::endl;

    BJReal ev = 0;
    for (int card = 1; card <= 10; ++card) {
        shoe.deal(card);
        BJPlayer* eor = new BJPlayer{shoe, rules, *basic, progress};
        shoe.reset();
        ev += eor->getValue() * shoe.getProbability(card);
    }
    std::cout << ev << std::endl;
}

with the following results (as fraction of initial wager):

-0.0033026803733750346
-0.0033026803733748255

These two values are close, but not exactly equal. Granted, these calculations were done using the limited precision of 64-bit floating-point arithmetic, with the accompanying rounding errors. But how can we tell whether the differences are due to limited numeric precision, or due to errors in the algorithm itself?

To answer this question, I recently updated the code to typedef the numeric data type (BJReal in the above example code) used for computing probabilities and expected values. This type defaults to double, but we can use the arbitrary-precision math::Rational instead– but otherwise using exactly the same algorithm– and verify that the expected return in the above scenario is exactly:

and that we get exactly this same result whether played from the full shoe or averaged over each possible burn card. In other words, the discrepancy in the lowest order digits above was due entirely to double-precision rounding error.

Of course, as shown in the high card “ties push” excursion described above, this doesn’t prove that the algorithm is correct. But inequality would have proved that the algorithm is incorrect. For this reason, I think this is a useful additional test to consider when evaluating algorithms like these.

Introduction

Alice and Bob are at it again, this time playing a game called Binary Search. Alice begins by secretly selecting an integer key from 1 to , then Bob makes a series of integer guesses, to each of which Alice responds indicating whether the guess is correct, or higher than her selected key value, or lower. Bob must pay Alice one dollar for each guess (including the final correct guess). How much should Alice be willing to pay to play this game with Bob?

This post is motivated by a recent blog post by John Graham-Cumming, discussed on Hacker News, about a similar game apparently used by Steve Ballmer in Microsoft job interviews. In the linked clip, Ballmer plays the role of Alice, and offers $6 to play the game. You’re Bob interviewing for a job: should you take the bet? Ballmer suggests that the answer is no, since as he puts it, “I [Alice/Steve] can pick numbers specifically that are hard for you [Bob] to get.”

As many have noted in the subsequent discussion, this comment certainly seems to allow for Alice/Steve picking numbers adversarially, that is, not necessarily uniformly at random. But presumably Bob is then allowed similar leeway, not being forced to execute textbook binary search that always guesses the midpoint of the interval of remaining possible numbers. So who has the greater strategic advantage? I’m not sure. I think this is a hard problem, as I’ll try to make the case below.

Solving smaller problems

I wrote some Mathematica code (on GitHub) to compute exact optimal strategies for Alice and Bob for smaller versions of this problem, up to . Alice’s pure strategies are simply the set of possible secret keys. Bob is more interesting; we can identify each possible pure search strategy as an in-order binary tree with vertices. The following function enumerates these (it’s always a nice surprise when my favorite sequence turns up in another unexpected place):

(* List all possible search strategies with key space[n]. *)
strategies[n_] := strategies[n] = Module[
   {guess},
   If[n === 0, {{}},
    Join @@ Table[Join @@ Outer[
        {guess, #1, #2} &,
        strategies[guess - 1],
        guess + strategies[n - guess],
        1
        ],
      {guess, 1, n}
      ]]]

Armed with Alice’s and Bob’s pure strategies, we can compute the matrix of payoffs resulting from each possible strategy match-up, then solve the corresponding linear programming problem to compute the expected value of the game and the Nash equilibrium mixed strategies for each player.

Results and conjectures

The results for up to 13 exhibit a consistent pattern suggesting that an optimal strategy for Alice is to pick numbers “almost uniformly” at random… but twice as likely to pick a number on the edge. More precisely, pick 1 or each with probability , and pick any other number with probability . I think it is likely that this is an optimal strategy for Alice in the case as well. (However, it’s interesting to observe that the presumed “Ballmer strategy” of picking uniformly from the subset of numbers for which the normal binary search would yield the maximum number of guesses has the same expected payoff.)

Bob’s strategy, on the other hand, is harder to describe in simple terms. Let’s consider for a specific example. Of the 742,900 possible pure search strategies, an optimal mixed strategy for Bob is to select from just the following baker’s dozen with the indicated probabilities, where each search strategy is specified as the pre-order traversal of the corresponding binary tree:

Prob.      Strategy
------------------------------------------------
0.381579   7  3  1  2  5  4  6 11  9  8 10 13 12
0.131579   8  4  2  1  3  6  5  7 12 10  9 11 13
0.131579   6  2  1  4  3  5 10  8  7  9 12 11 13
0.0625     9  4  1  2  3  6  5  8  7 12 10 11 13
0.060855   5  2  1  4  3 10  8  6  7  9 13 12 11
0.055921   9  4  2  1  3  6  5  8  7 12 10 11 13
0.055921   8  3  1  2  5  4  6  7 12 10  9 11 13
0.055921   5  2  1  4  3 10  8  6  7  9 12 11 13
0.055099   6  2  1  4  3  5 11  9  8  7 10 13 12
0.003289   9  5  2  1  4  3  7  6  8 12 10 11 13
0.003289   5  2  1  4  3  9  7  6  8 12 10 11 13
0.001645   5  2  1  4  3 10  7  6  9  8 13 12 11
0.000822   7  2  1  4  3  5  6 11  9  8 10 13 12

In other words, roughly 38% of the time Bob should conduct the usual even-split binary search, guessing the midpoint 7 to start. However, sometimes Bob might start with an initial guess as low as 5 or as high as 9.

In the original game where , I think it’s certainly the case that Bob should similarly be willing to start with guesses north of 51 or south of 50, but it’s unclear how extreme those departures might be.

Finally, what is Alice/Steve’s expected return in this game? Let’s set things up similarly to the Ballmer interview version of the problem, and suppose that Alice pays dollars to play, and wins back a dollar for each of Bob’s guesses. Then the following figure shows Alice/Steve’s expected return as a function of the number of possible keys to pick from.

My suspicion is that this sawtooth behavior persists, continuing to straddle zero expected return, so that it seems difficult to predict whether the particular game at has Alice in the red or the black, when both she and Bob use truly optimal strategies. This simple brute force computational approach certainly won’t work: the matrix of payoffs has 100 rows and roughly 10^57 columns.

Edit 2024-09-09: I was referred to Konstantin Gukov‘s very interesting article about this game, where they compute an explicit strategy for the case that proves that Alice/Ballmer’s expected return is indeed negative in that case. They do this by restricting Bob the Guesser’s set of available binary search strategies, from all possible strategies to just 586 carefully chosen ones, so that it is feasible to solve the corresponding linear programming problem. Since Alice’s expected return is negative against this handcuffed version of Bob, it must also be against a truly optimal guesser.

Gukov’s selection of subset of strategies is indeed well-chosen. I ran Gukov’s Python code to generate their search strategies for all up to 100, then computed the resulting expected return in each case. (Recall that we are assuming that Alice pays dollars to play, so that we can compare performance across a range of key space sizes.) The results are shown in the figure below, with Ballmer’s return against Gukov’s restricted guesser shown in red, and the performance against an optimal guesser shown in blue (essentially a reproduction of the earlier figure above).

This gives more evidence for our earlier speculation about how the advantage varies with , with a clearer picture of the pattern: it seems that Ballmer would have been safest by offering to play the game choosing numbers from 1 to 127– or 1 to 63, or 1 to 31, etc.– instead of 1 to 100. Powers of two are the worst for Ballmer, but as the key space increases from there, it eventually becomes advantageous again… until reaching the next power of two.

Years ago, I wrote a Python module using the VPython 3D graphics library to simulate a three-wheeled robot modeled after the Parallax Scribbler, that could execute Logo-like “turtle” motion commands. It worked really well as an educational tool for programming students, since using VPython’s graphics meant that you could manipulate the robot even from the interpreter prompt, a line at a time, and rotate and zoom the 3D view even while the robot is moving. No window event loop, no multithreading to think about.

I soon added a stall sensor and proximity sensors for detecting and navigating around obstacles, such as the walls of a maze as shown in the screenshot below.

from vturtle import *
robot = Robot(obstacles=maze())

In the intervening 14-ish years, VPython has undergone a significant refactoring. I’ve recently updated the robot simulator code accordingly, so that it now works in the currently latest VPython version 7.6.5. The code is on GitHub.