The 538 Riddler:New Casino Game

Welcome to my weekly look at 538’s The Riddler. Happily my answer to last week’s riddle proved to be correct once again, and I was happy to see some commenters who came to the same conclusions. Thanks for reading and commenting.

Due to a lack of time with Easter and all, I’m not very confident of my answer to this week’s Riddler, but I’ll let you decide. Here’s the riddle:

Suppose a casino invents a new game that you must pay $250 to play. The game works like this: The casino draws random numbers between 0 and 1, from a uniform distribution. It adds them together until their sum is greater than 1, at which time it stops drawing new numbers. You get a payout of $100 each time a new number is drawn.

For example, suppose the casino draws 0.4 and then 0.7. Since the sum is greater than 1, it will stop after these two draws, and you receive $200. If instead it draws 0.2, 0.3, 0.3, and then 0.6, it will stop after the fourth draw and you will receive $400. Given the $250 entrance fee, should you play the game?

Specifically, what is the expected value of your winnings?

Okay, when I read that part it seemed like this was going to be a no-brainer, but then Ollie added:

[Editor’s note: You know how sometimes your boss/teacher asks you not to bring your laptop to a meeting/class, and you’re all like, “Jeez, how am I gonna tweet during the meeting/class, then?” Well, I’m going to pull rank and ask that you don’t bring laptops to this Riddler. Next week I will highlight — nay, celebrate! — the elegant, pencil-and-paper solutions over the brutal, cold, silicon ones. We’re all on the honor system with one another, right? RIGHT?]

No computers! Pencil and paper? Egads. Okay, so I did a little bit of working this out and I came up with an answer I thought might be right. First of all, I assumed that you are guaranteed to see at least 2 numbers, because with all the real numbers between 0 and 1, the chances of getting 1 exactly are one over infinity, or zero. So the minimum amount you are going to lose in this game is $50. And it seems to me that the chances of seeing a third number are 50%, since the average value of the first number will be 0.5 and the chances that the second number will be greater than 0.5 is… 0.5, or 50%.

Okay, so right away we see that there’s a 50% chance we lose $50 (-$250 + $200), and a 50% chance that we win $50 (-$250 + $300), with the possibility of more money to come. This is a wonderful casino game; I would very much like to find the casino that offers it. A game like this might explain how somebody could bankrupt a casino. (Sad!)

So we have a 50% chance of seeing only 2 numbers and losing $50. What are the chances of seeing only 3 numbers? Well, I figure if 2 random numbers are below 1, then there’s a 2/3 chance that the 3rd number will break the ceiling. So 2/3 of what is remaining (1/2) is 1/3. So 1/2 the time we see 2 and only 2 numbers and 1/3 of the time we see 3 and only 3 numbers, and now we have 1/6 (1/2 – 1/3) remaining.

Chance of only 2 numbers: 1/2. 1 – 1/2 = 1/2 remaining.
* 2/3  =
Chance of only 3 numbers: 1/3. 1/2 – 1/3 = 1/6 remaining.
* 3/4 =
Chance of only 4 numbers: 1/6 * 3/4 = 1/8. 1/6 – 1/8 = 1/24 remaining.
* 4/5 =
Chance of only 5 numbers: 1/24 * 4/5 = 1/30. 1/24 – 1/30 = 1/120 remaining.
* 5/6 =
Chance of only 6 numbers: 1/120 * 5/6 = 1/144.

And so on.

So to get the Expected Value (EV) we can take the odds of each payout * its value and get:

$-50/2 + $50/3 + $150/8 + $250/30 + $350/144 = $21.18

Plus some increasingly smaller numbers for seeing 7, 8, 9, 10? numbers.

Ok, now that we got our EV answer, let’s cheat! (If you don’t want to see whether a computer confirms or disproves my EV answer, stop reading here.)

I wrote some code in Go (are you sick of Go yet, because I am not.)

package main

import (
	"fmt"
	"math/rand"
)

func doFloatTrial() int {
	total := 0.0
	numbers := 0
	
	for {
		numbers++
		total += rand.Float64()
		if total >= 1.0 {
			break;
		}
	}
	return numbers
}


func main() {

	trials := 1000000
	
	var results []int
	
	cash := 0.0
	
	for i := 0; i < trials; i++ {
		cash -= 250.0
		numbers := doFloatTrial()
		cash += 100.0 * float64(numbers)
		for a := len(results); a <= numbers; a++ {
			results = append(results, 0)
		}	
		results[numbers]++
	}

	fmt.Printf("EV = %.3f\n", cash / float64(trials));
	for i := 0; i < len(results); i++ {
		fmt.Printf("Results[%d] = %d (%.2f) \n", i, results[i], float64(trials)/float64(results[i]))
	}

}

 

And here’s the output:

EV = 21.826
Results[0] = 0 (+Inf) 
Results[1] = 0 (+Inf) 
Results[2] = 499604 (2.00) 
Results[3] = 334041 (2.99) 
Results[4] = 124717 (8.02) 
Results[5] = 33341 (29.99) 
Results[6] = 6939 (144.11) 
Results[7] = 1157 (864.30) 
Results[8] = 184 (5434.78) 
Results[9] = 15 (66666.67) 
Results[10] = 2 (500000.00)

Hurray! Our EV from the code is $21.83 and we can see that the odds of drawing each count of numbers is (with variance) what we had calculated by hand. Math! Follow this link to run the code yourself in your browser.

Thanks for reading. I’m not 100% sure of this one so I look forward to better approaches to this problem in the comments.

The 538 Riddler:Weird Guy in Trench Coat

Hello, and welcome to my look at 538’s Riddler. Last week I managed to get the correct answer using a Monte Carlo, and Go’s really nice NormFloat64 function. Commenters squareandrare, mathisallaround, and guyshavit all shared voodoo mathy ways to get the answer, which I appreciated, including the chi-squared distribution. Thanks for reading and commenting!

This week’s Riddler concerns talking to strange men in trench coats.

A man in a trench coat approaches you and pulls an envelope from his pocket. He tells you that it contains a sum of money in bills, anywhere from $1 up to $1,000. He says that if you can guess the exact amount, you can keep the money. After each of your guesses he will tell you if your guess is too high, or too low. But! You only get nine tries. What should your first guess be to maximize your expected winnings?

Obviously if a man in a trench coat offers you money the answer is to run. But let’s answer the question anyway. My immediate first guess for this was $744, which was wrong. But only because I’m pretty sure the answer is $745.

How did I get this number? Well, the guessing game where someone tells you whether your guess is higher or lower is a classic example of a binary search, which is one of the most important algorithms in computer science. With 9 guesses, we know that we can guarantee coverage of 511 (2^9 – 1) values. That leaves 489 values for which we are going to not be able to guess the right answer. Sad!

Now the question is, which values do we want to cover? Well, without any knowledge to the contrary we should assume that all values between $1 and $1000 are equally likely, assuming our trench-coated friend is trustworthy (and aren’t they all?). So since there is just as much chance that there is $1000 in the envelope as $1, we might as well make sure we cover the top 511 values.

Put another way, imagine the man in the trench coat offers you only one guess. What should you do? (Run.) Well, you have just as much chance of being right by guessing $1000 as $1, so guess $1000 and hope for the best. Your expected value (EV) is the odds of being right times the payout, so if you guess $1000, your EV is $1. If you guess $1, your EV is $.001, which is… less.

So covering the top 511 values would be everything from $490 to $1000. The middle value between those two numbers is $745. I wrote some code to show this, and then ran 1000 simulations for each amount of money that is in the envelope. Here it is in Go:

package main
// Jon Wiesman, somedisagree.com

import "fmt"

func binSearchRecursive(min, max, value, guesses int) bool {

	if guesses == 0 {
		return false
	}
	if min > max {
		return false
	}

	guess := (min + max) / 2
	if guess == value {
		return true
	}

	if guess < value {
		return binSearchRecursive(guess+1, max, value, guesses-1)
	} else {
		return binSearchRecursive(min, guess-1, value, guesses-1)
	}

	return false
}

func binSearchIterative(min, max, value, guesses int) bool {
	for ; guesses > 0; guesses-- {
		if min > max {
			return false
		}

		guess := (min + max) / 2
		if guess == value {
			return true
		}

		if guess < value {
			min = guess + 1
		} else {
			max = guess - 1
		}

	}
	return false
}

func DoTrialRecursively(min, max int) {
	trials := 1000
	found := 0
	totalEarned := 0

	for i := 1; i <= 1000; i++ {
		if binSearchRecursive(min, max, i, 9) {
			found++
			totalEarned += i
		}
	}

	fmt.Printf("Recursively, with first value = %d, %d/%d success, total earned = %d\n", (min+max)/2, found, trials, totalEarned)
}

func DoTrialIteratively(min, max int) {
	trials := 1000
	found := 0
	totalEarned := 0

	for i := 1; i <= 1000; i++ {
		if binSearchIterative(min, max, i, 9) {
			found++
			totalEarned += i
		}
	}

	fmt.Printf("Iteratively, with first value = %d, %d/%d success, total earned = %d\n", (min+max)/2, found, trials, totalEarned)
}

func main() {
	DoTrialIteratively(489, 1000)
	DoTrialIteratively(490, 1000)
	DoTrialRecursively(490, 1000)
	DoTrialRecursively(491, 1000)

}

 

And here’s a link to the Go playground for this code where you can play with the numbers and run it in your browser. Doing so gives this output:

Iteratively, with first value = 744, 511/1000 success, total earned = 380184
Iteratively, with first value = 745, 511/1000 success, total earned = 380695
Recursively, with first value = 745, 511/1000 success, total earned = 380695
Recursively, with first value = 745, 510/1000 success, total earned = 380205

So your EV if you start with a first guess of $745 is $380.695, and you’ve got a 51.1% chance of making some money. Pretty good deal. (You should still run.)

My binary search function differs from most binary searches in one small way. Here, I pass in the number of guesses and decrement it with each recursive call. When it reaches 0, I return false. I used recursion here first because it’s easiest for me to write but iteration might be easier for non CS people to follow so I wrote it iteratively as well. Both work.

If you call either function with min and max at $1 and $1000, respectively, giving a first guess of $500, you will still succeed 511 times but your expected value would only be $255.175. This is because the 489 missed values would be evenly distributed instead of all at the bottom. Try it.

So, that’s my answer: $745. I believe it is right unless there is some method to guarantee more than 511 found values with only 9 guesses. As always, I’m always willing to admit I am wrong so if I am, I’d love to know.

Thanks for reading.

 

The 538 Riddler: Free Throw Grannies

Hello, and welcome to this week’s look at the 538 Riddler. Last week concerned optimal strategy for a (very) simple game show, and happily my answer was correct. This week’s Riddler concerns free throws again.

Consider the following simplified model of free throws. Imagine the rim to be a circle (which we’ll call C) that has a radius of 1, and is centered at the origin (the point (0,0)). Let V be a random point in the plane, with coordinates X and Y, and where X and Y are independent normal random variables, with means equal to zero and each having equal variance — think of this as the point where your free throw winds up, in the rim’s plane. If V is in the circle, your shot goes in. Finally, suppose that the variance is chosen such that the probability that V is in C is exactly 75 percent (roughly the NBA free-throw average).

But suppose you switch it up, and go granny-style, which in this universe eliminates any possible left-right error in your free throws. What’s the probability you make your shot now? (Put another way, calculate the probability that |Y| < 1.)

So I’ve been doing all these Riddler exercises in C++ but this week I thought I’d mix it up and give the Go language a shot. I’ve never written anything in Go but a friend of mine sent me this text a while back:

text_from_sean

As you can see, I have an exciting social life. So anyway I’ve been wanting to try it out for a while and I figured this blog would be a good opportunity. I’ve been doing the Go tutorials that you can find here.

Another benefit of doing the code (when possible) in Go is that I can link to a version of my code that you can run yourself in the browser, and then make changes to see the impact of those changes. I’m not sure that anyone was copying my C++ code into a compiler, creating a project, building, then running it. But this way running the code is literally just a click away.

Okay, so I thought this week’s Riddler was pretty straightforward from a “getting an answer” point of view. (Which probably means I’m totally wrong.) If I simulate a normalized random distribution of 2D points with some variance such that 75% of those points fall inside a circle with radius 1, then what would be the probability of falling in the circle if the x-component was always 0?

Luckily, Go has a really nice NormFloat64 function that generates normalized random numbers with a standard deviation of 1. With such variance, about 40% of random points will be inside the circle, using the good old Pythagorean Theorem. If I understand the problem correctly all I need to do is find the standard deviation that will result in 75% of points falling inside the circle. Then, I can calculate the magnitude of points where the x-component is 0, and find how frequently they fall in the circle.

So the big task here was to find the right value for the deviation factor. I did this with trial and error and came up with 0.60054. I’m sure there’s a math way to find this value (or a better value) and I look forward to seeing it (Attn: Alex!). Once you start generating x- and y-components, find their magnitude using Sqrt(x*x + y*y) and see if it is below 1. If it is, then that’s a made shot. Also, check to see what would happen if the x-component is 0. You can do that as Sqrt(0*0 + y*y) or Sqrt(y*y) or Abs(y).

Here’s the code in Go:

// 538 Riddler: Should You Shoot Free Throws Underhand?
// http://fivethirtyeight.com/features/should-you-shoot-free-throws-underhand/
// Jon Wiesman, somedisagree.com
package main

import (
	"fmt"
	"math"
	"math/rand"
)

func main() {

	overhand := 0
	granny := 0

	trials := 1000000
	r := rand.New(rand.NewSource(int64(538)))

	dev := 0.60054
	
	for i := 0; i < trials; i++ {
		x := r.NormFloat64() * dev
		y := r.NormFloat64() * dev

		d := math.Sqrt(x*x + y*y)
		g := math.Abs(y)

		if d < 1 {
			overhand++
		}
		if g < 1 {
			granny++
		}
	}

	fmt.Printf("Overhand = %f%%\n", float64(overhand)*100.0/float64(trials))
	fmt.Printf("Granny = %f%%\n", float64(granny)*100.0/float64(trials))
}

 

And here’s a link that will take you to a Go playground where you can run that code and see these results:

Overhand = 75.000800%
Granny = 90.439500%

Try it! You can also play around with different random seed values (I used 538, natch) and see what happens. I think the actual make percentage for Granny-style is more like 90.41%, depending on random seed and number of trials.

Thanks for reading!

The 538 Riddler: Encryption Blues

This week’s Riddler was mostly fun with a solid day of frustration thrown in. I was able to get the answers but I didn’t use much code to do it, but I’ll show what I did write. Also, a brief note about scheduling for the future. I’m going to start posting my write-ups after the 9 PM (PST) deadline for submitting answers on the Riddler site. I’ll still submit my answers to them on time but I don’t want to post answers here before the deadline. There’s a five-day window between the submission deadline and the next riddle and I think it makes sense to post answers in that window so as not to spoil anyone who might be looking for hints, but not interested in finding the answer. Once the submission deadline has passed, I’ll go ahead and post my answer.

Here’s this week’s problem:

Here are four different coded messages. Each has a different encryption scheme and they progress, I think, from easiest to toughest to crack. Submit the decrypted messages as your answers.

1. A zsnw kmuuwkkxmddq kgdnwv lzw XanwLzajlqWayzl Javvdwj!

2. xckik acvlbeg oz mmqn xnlautw. gzag, mwcht, kbjzh… ulw cpeq edr mom dhqx lksxlioil?

3. hy vg nw rh ev pr is or tf?

4. 😎😊, 😓😇😀😓’😒 😈😓. 😍😎😖 😆😄😓 😁😀😂😊 😓😎 😖😎😑😊. 😇😄😘, 😀😓 😋😄😀😒😓 😈😓’😒 😅😑😈😃😀😘.

So the first two codes I was able to get in very little time. I suspected that the first one would be a simple Caesar, (aka Wheel or Rotation) cipher. I wrote a simple function that decrypts characters by shifting letter values based on a shift parameter, and wrapping if necessary.

#define ABC_SIZE 26

wchar_t WheelDecodeChar(wchar_t ch, int shift)
{
    wchar_t decoded = 0;

    if(ch >= 'a' && ch <= 'z')
    {
        decoded = ch + shift;
        while(decoded > 'z')
            decoded -= ABC_SIZE;
        while(decoded < 'a')
            decoded += ABC_SIZE;
    }
    else if(ch >= 'A' && ch <= 'Z')
    {
        decoded = ch + shift;
        while(decoded > 'Z')
            decoded -= ABC_SIZE;
        while(decoded < 'A')
            decoded += ABC_SIZE;
    }
    else
    {
        decoded = ch;
    }
    return decoded;
}

void WheelDecode(const wchar_t *message, int shift)
{
    wchar_t decoded[256] = {};
    int len = wcslen(message);
    for(int ch = 0; ch < len; ch++)
    {
        decoded[ch] = WheelDecodeChar(message[ch], shift);
    }
    wprintf(L"Possible message = %s\n", decoded);
}

void WheelDecodeAll(const wchar_t *message)
{
    for(int i = 0; i < ABC_SIZE; i++)
    {
        WheelDecode(message, i);
    }
}

 

When I called WheelDecodeAll I got this output:

Possible message = A zsnw kmuuwkkxmddq kgdnwv lzw XanwLzajlqWayzl Javvdwj!
Possible message = B atox lnvvxllyneer lheoxw max YboxMabkmrXbzam Kbwwexk!
Possible message = C bupy mowwymmzoffs mifpyx nby ZcpyNbclnsYcabn Lcxxfyl!
Possible message = D cvqz npxxznnapggt njgqzy ocz AdqzOcdmotZdbco Mdyygzm!
Possible message = E dwra oqyyaoobqhhu okhraz pda BeraPdenpuAecdp Nezzhan!
Possible message = F exsb przzbppcriiv plisba qeb CfsbQefoqvBfdeq Ofaaibo!
Possible message = G fytc qsaacqqdsjjw qmjtcb rfc DgtcRfgprwCgefr Pgbbjcp!
Possible message = H gzud rtbbdrretkkx rnkudc sgd EhudSghqsxDhfgs Qhcckdq!
Possible message = I have successfully solved the FiveThirtyEight Riddler!
Possible message = J ibwf tvddfttgvmmz tpmwfe uif GjwfUijsuzFjhiu Sjeemfs!
Possible message = K jcxg uweeguuhwnna uqnxgf vjg HkxgVjktvaGkijv Tkffngt!
Possible message = L kdyh vxffhvvixoob vroyhg wkh IlyhWkluwbHljkw Ulggohu!
Possible message = M lezi wyggiwwjyppc wspzih xli JmziXlmvxcImklx Vmhhpiv!
Possible message = N mfaj xzhhjxxkzqqd xtqaji ymj KnajYmnwydJnlmy Wniiqjw!
Possible message = O ngbk yaiikyylarre yurbkj znk LobkZnoxzeKomnz Xojjrkx!
Possible message = P ohcl zbjjlzzmbssf zvsclk aol MpclAopyafLpnoa Ypkksly!
Possible message = Q pidm ackkmaancttg awtdml bpm NqdmBpqzbgMqopb Zqlltmz!
Possible message = R qjen bdllnbboduuh bxuenm cqn OrenCqrachNrpqc Armmuna!
Possible message = S rkfo cemmoccpevvi cyvfon dro PsfoDrsbdiOsqrd Bsnnvob!
Possible message = T slgp dfnnpddqfwwj dzwgpo esp QtgpEstcejPtrse Ctoowpc!
Possible message = U tmhq egooqeergxxk eaxhqp ftq RuhqFtudfkQustf Duppxqd!
Possible message = V unir fhpprffshyyl fbyirq gur SvirGuveglRvtug Evqqyre!
Possible message = W vojs giqqsggtizzm gczjsr hvs TwjsHvwfhmSwuvh Fwrrzsf!
Possible message = X wpkt hjrrthhujaan hdakts iwt UxktIwxginTxvwi Gxssatg!
Possible message = Y xqlu ikssuiivkbbo ieblut jxu VyluJxyhjoUywxj Hyttbuh!
Possible message = Z yrmv jlttvjjwlccp jfcmvu kyv WzmvKyzikpVzxyk Izuucvi!

So, we can see the 8th line is the answer we are looking for.

The second question is a little harder, but I was able to get it pretty easily. From looking at the string, “xckik acvlbeg oz mmqn xnlautw. gzag, mwcht, kbjzh… ulw cpeq edr mom dhqx lksxlioil?” you can kind of tell that the encryption alphabet must be changing. I was able to guess that this was probably a Vigenère cipher and wrote a function to modify the shift value passed into WheelDecodeChar based on a key. This works by using the first letter of the key as the shift value for the first letter, then the second letter in the key for the second letter, and so on. If the key is shorter than the message, you just go back to the beginning of the key.

void WheelDecodeKey(const wchar_t *message, const wchar_t *key)
{
    int keyLen = wcslen(key);
    int msgLen = wcslen(message);

    wchar_t decoded[256] = {};
    int keyIdx = 0;
    for(int i = 0; i < msgLen; i++)
    {
        int shift = tolower(key[keyIdx % keyLen]) - 'a';
        decoded[i] = WheelDecodeChar(message[i], -shift);
        if(IsLetter(message[i]))
            keyIdx++;
    }
    wprintf(L"Possible message = %s\n", decoded);
}

 

The only trick was finding which key to use. I tried “Riddler”, “TheRiddler”, and then finally tried “FiveThirtyEight” and got this:

Possible message = super tuesday is this tuesday. cruz, trump, rubio... who will
 win the most delegates?

That looks about right. Cool, so I got the first two in pretty short order. How hard could the third one be?

What followed was about a full day of frustration before I decided to move on to the fourth one. The fourth one uses emojis instead of letters but I decided to fix that by making a function to assign letters to the emojis.

void MakeLettersFromEmojis(const wchar_t *message, wchar_t *out)
{
    // SUPER UGLY CODE with lots of hard-coded assumptions

    map<unsigned int, wchar_t> mapping;

    int len = wcslen(message);
    wchar_t ch = 'a';
    bool esc = false;
    int outIndex = 0;

    for(int i = 0; i < len; i++)
    {
        wchar_t m = message[i];
        if(m == ' ' || m == '.' || m == ',' || m == '\'' || m == 0x2019)
        {
            if(m == 0x2019)
                m = '\'';
            out[outIndex] = m;
            outIndex++;
            continue;
        }

        unsigned int m2 = m;
        if(esc)
        {
            m2 = (u16esc << 16) | m;
            esc = false;
        }
        else if(m == u16esc)
        {
            esc = true;
            continue;
        }

        auto itr = mapping.find(m2);
        if(itr == mapping.end())
        {
            mapping[m2] = ch;
            out[outIndex] = ch;
            outIndex++;
            ch++;
        }
        else
        {
            out[outIndex] = itr->second;
            outIndex++;
        }
    }
    out[outIndex] = 0;
    wprintf(L"%s\n", out);
}

 

That produced the following output:

ab, cdec'f gc. hai jkc lemb ca ianb. dko, ec pkefc gc'f qngreo.

Now this isn’t a simple Caesar cipher. This is a completely random alphabet. There are some really interesting algorithms that use frequency analysis to come up with likely cipher alphabets, but I just solved this by hand. (This looks like a lot of fun to code; I might revisit it later.) With the apostrophes, I guessed that the “f” was an “s”, which led me to guess that the “c” was a “t”.

ab, tdet's gt. hai jkt lemb ta ianb. dko, et pkest gt's qngreo.

That led me to guess that "de" in the second word was "ha" to spell "that's". Also the "g" would be an "i". This leads to:
ab, that's it. hai jkt lamb ta ianb. hko, at pkast it's qnirao.

So all those guesses were correct, which is not always how solving cryptograms go, but I got lucky. From there it was a matter of some trial and error to finally get:

ok, that's it. now get back to work. hey, at least it's friday.

Okay, now back to the third message and my frustration.

I spent the better part of a day trying to figure out what the encryption was for this string:

hy vg nw rh ev pr is or tf?

I decided it had to be digraphs and there’s an interesting cipher called a digraph substitution cipher that I coded up. I wrote code to generate all 676 possible digraph substitution alphabets based on shifting rows and columns by 26 characters. Then I used the output of those and scored each possible string by adding up the resultant digraphs and their frequency from this table. That was really cool except that this string wasn’t encrypted with a digraph substitution cipher, so none of the 676 possible alphabets produced any good candidates.

Finally I was chatting with a friend about it who asked, “have you tried Playfair?” The answer was no because I had never heard of Playfair. Here’s the wikipedia link for the Playfair cipher. It’s pretty cool, and JFK used it once!

Anyway, after trying a plain alphabet I used the key, “FIVETHRYG” (no repeated letters) and got the following output:

ar ey ou ha vi ng fu ny et?

(are you having fun yet?)

And that’s it! If I had had more time I probably would have coded that up, which would have been fun, but I spent too much time on the digraphs.

Thanks for reading.

The 538 Riddler: Air Travel Nightmare

So last week’s Riddler was a good news/bad news situation. The bad news is that I was wrong, although I did illustrate the correct answer in a follow-up post. The good news is that I still managed to get a shout-out from Ollie in this week’s Riddler, even earning the honorary title of FotR, which I hope means Friend of the Riddler.

Anyway, this week’s Riddler concerns the worst person alive and a full airplane.

There’s an airplane with 100 seats, and there are 100 ticketed passengers each with an assigned seat. They line up to board in some random order. However, the first person to board is the worst person alive, and just sits in a random seat, without even looking at his boarding pass. Each subsequent passenger sits in his or her own assigned seat if it’s empty, but sits in a random open seat if the assigned seat is occupied. What is the probability that you, the hundredth passenger to board, finds your seat unoccupied?

So the answer is 1/2, half, 50%, however you want to say it. And this is true for the last passenger no matter how many passengers are on a (full) plane. I actually worked this out before writing any code (but don’t worry, I still wrote the code!) It’s trivial to show for the case where there are only two passengers: there are two seats, and the person in front of you is the worst person alive. He either takes his seat, or yours, randomly. Fifty percent, done.

For three people (1. Worst Person Alive, 2. Rando Calrissian, and 3. You) and three seats, we can still do it in our heads. WPA takes one of three seats, so right away you have a 1/3 chance of not getting your seat. But there is also a 1/3 chance that he takes Rando’s seat, in which case Rando will randomly take one of the two remaining seats, yours or WPA’s assigned seat. One-third times one-half is 1/6, so there is a 1/6 chance you will find Rando sitting in your seat, and a 1/3 chance WPA will be sitting there. So 1/3 + 1/6 = 1/2. Boom, done.

Once I saw that the probability remained at 50% for three passengers, I suspected that it would be 50% for any number of passengers. I wrote a Monte Carlo to show that.

// AirplaneSeating.cpp : Defines the entry point for the console application.
// www.somedisagree.com, Jon Wiesman

#include "stdafx.h"
#include <vector>

#define SEAT_COUNT 100
#define EVIL_COUNT 1
#define AVAILABLE -1

class Plane
{
    std::vector<int>    m_seats;
    int                 m_available;

public:
    Plane(int seats)
    {
        m_seats.resize(seats);
        for(int i = 0; i < seats; i++)
            m_seats[i] = AVAILABLE;
        m_available = seats;
    }

    void    ClaimSeat(int seat, int passenger) 
    {
        m_seats[seat] = passenger;
        m_available--;
    }

    int     FindUnclaimedSeat() const
    {
        if(m_available == 0)
            return -1;
        int nth = rand() % m_available;

        int empty = 0;
        for(size_t i = 0; i < m_seats.size(); i++)
        {
            if(IsSeatAvailable(i))
            {
                empty++;
                if(empty > nth)
                    return i;
            }
        }
        return -1;
    }

    bool    IsSeatAvailable(size_t seat) const
    {
        return (seat >= 0 && seat < m_seats.size() && m_seats[seat] == AVAILABLE);
    }

    void    GetSeat(int passenger, bool evil)
    {
        if(evil || !IsSeatAvailable(passenger))
        {
            int seat = FindUnclaimedSeat();
            ClaimSeat(seat, passenger);
        }
        else
        {
            ClaimSeat(passenger, passenger);
        }
    }

    void    SimulateBoarding(size_t evilCount)
    {
        for(size_t i = 0; i < m_seats.size(); i++)
        {
            GetSeat(i, (i < evilCount));
        }
    }

    void    CompileStats(std::vector<int> &stats, int &wrongSeats)
    {
        for(size_t i = 0; i < m_seats.size(); i++)
        {
            if(m_seats[i] == i)
                stats[i]++;
            else
                wrongSeats++;
        }
    }
};

int _tmain(int argc, _TCHAR* argv[])
{
    std::vector<int> stats;
    int wrongSeats = 0;
    stats.resize(SEAT_COUNT);

    int trials = 100000;

    for(int i = 0; i < trials; i++)
    {
        Plane plane(SEAT_COUNT);

        plane.SimulateBoarding(EVIL_COUNT);
        plane.CompileStats(stats, wrongSeats);
    }

    printf("After %d trials with %d evil passengers:\n", trials, EVIL_COUNT);
    for(int i = 0; i < SEAT_COUNT; i++)
    {
        printf("Passenger %d seated properly %.2f%% of time\n", i, 100.0 * (double)stats[i]/(double)trials);
    }
    printf("Average of %.2f wrongfully seated passengers\n", (double)wrongSeats / (double)trials);

	return 0;
}

 

This code (download) is pretty simple. The Plane class keeps an array of seats and who (which passenger) is seated in each seat. It also keeps the number of empty seats at any time to make it easier to quickly choose a random seat. (For the purpose of this simulation it was convenient to say that the passenger index and seat assignment are the same.)

The constructor initializes the seat array to all available and sets the available count to the seat count. SimulateBoarding tells each passenger to find a seat. Finding a seat is easy. If you’re not evil and your seat is available, sit down. Otherwise, find a random available seat. Either way, call ClaimSeat to sit down and decrease the available member variable. A CompileStats function will increment the number of correct passenger-seat assignments and keep a running count of wrong ones. Finally we run this simulation 100,000 times and print out the results.

After 100000 trials with 1 evil passengers:
Passenger 0 seated properly 1.00% of time
Passenger 1 seated properly 99.03% of time
Passenger 2 seated properly 99.00% of time
Passenger 3 seated properly 99.00% of time
Passenger 4 seated properly 99.01% of time
Passenger 5 seated properly 98.93% of time
Passenger 6 seated properly 98.97% of time
Passenger 7 seated properly 98.95% of time
Passenger 8 seated properly 98.92% of time
Passenger 9 seated properly 98.98% of time
Passenger 10 seated properly 98.87% of time
.....
Passenger 90 seated properly 90.86% of time
Passenger 91 seated properly 90.03% of time
Passenger 92 seated properly 88.83% of time
Passenger 93 seated properly 87.42% of time
Passenger 94 seated properly 85.39% of time
Passenger 95 seated properly 83.11% of time
Passenger 96 seated properly 79.98% of time
Passenger 97 seated properly 74.74% of time
Passenger 98 seated properly 66.45% of time
Passenger 99 seated properly 49.95% of time
Average of 5.19 wrongfully seated passengers

So I ran this code really just to confirm that my guess of 50% would be correct. However, I saw something which kind of jogged my memory. The average number of wrongfully seated passengers out of 100 was 5.19. This was the average number of traffic jams with 100 cars from two weeks ago! (The exclamation point at the end of the last sentence might be the nerdiest exclamation point ever.) It turns out that the number of wrongfully seated passengers follows a Harmonic Series.

What this means is that the odds of anyone getting a wrong seat follows this pattern, starting from the last passenger to board: 1/2 + 1/3 + 1/4 + … You can see from the output that the chances of each passenger not getting their seat is indeed following that pattern, with a little Monte Carlo variance thrown in.

“But wait,” you say, probably not, but let’s pretend you do, “a Harmonic Series starts with one and you don’t have a one in there.” Oh but I do. But it comes at the end. Remember that in our simulation, Worst Person Alive boards first and therefore has a 99/100 chance of sitting in the wrong seat. The next person to board, let’s call him Rando Paul, has a 1/100 chance of finding WPA in  his seat. 99/100 + 1/100 = 1. The next person then has a 1/99 chance and on down to 1/2 for the last passenger. So the average number of wrongly seated passengers for N passengers with the Worst Person Alive seating first boils down to 1 + 1/2 + 1/3 + … + 1 / (N – 1).

Traffic jams and airline travel: unexpected sources of harmony. Thanks for reading!

 

Follow up to Duck vs. Dog

Okay, well, my luck has officially run out. I was definitely, most certainly wrong wrong wrong about how much faster the dog could travel and still allow the duck to escape. Thanks once again to commenters Alex Mahdavi and Lego Haryanto who found the right answer. Alex pointed me to this DataGenetics blog post from a few years ago that answered the question with a monster and a rower in a boat. As I suspected the answer was more mathy than my simulation could handle, although a lot of it was algebra and trig, with a bit of calculus.

The nice thing was that once I understood the solution, it was really easy to rewrite the SmartDuckStrategy to use the “J” strategy outlined in the post. The result is exactly what was predicted by Datagenetics. Here’s a GIF:

duckvdog_Correct

Yay! It’s nice when things work out, even if I didn’t solve this one (or get close) myself. The new code is a little simplistic in that it assumes the dog will start at a position of 0 degrees, but here it is:

DPoint SmartDuckStrategy::GetDestinationVector(Pond *pond)
{
    if(pond->GetBeeline())
        return pond->GetLastDuckVector();

    DPoint focus(0, -1 / (pond->GetDogSpeed() * 2.0));

    DPoint duckPos = pond->GetDuckPos();

    DPoint delta = duckPos - focus;

    double angleFromFocus = atan2(delta.y, delta.x);

    if(angleFromFocus > -PI / 2.0 && angleFromFocus < 0)
    {
        pond->SetBeeline();
        return DPoint(1, 0);
    }
    angleFromFocus += PI / 2.0;

    DPoint vec(cos(angleFromFocus), sin(angleFromFocus));

    return vec;
}

So I take some comfort that my underlying simulation was solid enough that just changing the behavior of the duck was sufficient to illustrate the right answer. It’s not much, but it’s all I have!

Thanks again to Alex and Lego for posting comments. It means a lot.

The 538 Riddler: A Dog and a Duck in Purgatory

UPDATED 2/15 2:09 PM: Gonna save you the suspense on this one, folks. I was wrong, very wrong. See the follow up here.

Things are turning a bit dark at the 538 Riddler as this week Ollie imagines a frantic dog and a beleaguered duck trapped in a never-ending Kafkaesque chase. My answer for last week’s puzzle was correct, but as I said later, it was much more complicated than it needed to be.

This week, however, I don’t think I will even get a correct-but-overly-complex answer. I’m positive that my answer is wrong.

Here’s the scenario for these poor creatures:

There is a duck paddling in a perfectly circular pond, and a Nova Scotia duck tolling retriever prowling on the pond’s banks. The retriever very much wants to catch the duck. The dog can’t swim, and this particular duck can’t take flight from water. But the duck does want to fly away, which means it very much wants to make it safely to land without meeting the dog there. Assume the dog and the duck are both very smart and capable of acting perfectly in their best interests — they’re both expert strategists and tacticians. Say the duck starts in the center of the pond, and the dog somewhere right on the pond’s edge. To ensure that the duck could never escape safely, how many times faster would the dog have to be able to run relative to how fast the duck can swim? (Assume the dog is a good boy, yes he is, and thus very patient. The duck can’t simply wait him out. Also assume both animals can change direction without losing speed, and that the duck can take flight safely immediately after reaching shore if the dog isn’t there.)

Okay, so obviously there’s a mathy solution to this problem, but I don’t know it. I wrote a program to simulate the dog and duck and gave the duck some different strategies for getting to the shore. One thing that is obvious: the dog will never catch the duck. Either the duck will escape because the dog is too slow, or the dog will circle the pond forever as the duck swims back in forth in vain. Geez, Ollie, next to this, last week’s nightmare traffic problem is looking downright cheerful.

Designing the simulation was pretty fun, even if I don’t think it led me to the right solution. I basically came up with four different strategies that the duck might take to escape the pond. I call these Dumb, Cowardly, Clever, and Smart. Each strategy has a threshold at which it stops working, which is question being asked.

The Dumb strategy simply involves moving due west and hoping the dog is too slow. It will succeed as long as the dog is moving slower than PI * the duck’s speed.

duckvdog_Dumb

The Cowardly strategy just tries to avoid the dog at all costs, changing the duck’s direction to be away from the dog. Interestingly, the Cowardly strategy does worse than the dumb strategy for the duck! The dog can “scare” the duck into remaining in the pond by getting close enough to scare it back towards the center of the pond, even when moving slower (0.9483) than PI * the duck’s speed. All ducks die, but a cowardly duck is doomed to swim in a pond forever or something.

duckvdog_Cowardly

The Clever strategy works like the Cowardly strategy except that the duck is constantly calculating how long it will take to reach the shore in a straight line, and whether the dog can make it there first. If the duck can make it to the edge first, it makes a beeline for the shore. The Clever strategy requires the dog to travel 1.2093 * PI * the duck’s speed to prevent takeoff.

duckvdog_Clever

Finally, my last strategy was the Smart strategy. This is what I submitted for an answer, but I will be shocked if it is right. The Smart strategy acts like the Clever strategy except the duck will scan a range of spots on the shore and if ANY of them can be reached before the dog, it will beeline. It represents only a small improvement over Clever, requiring the dog to travel 1.2552 * PI * the duck’s speed to trap the fowl.

duckvdog_Smart

A couple of things to note here. When the duck escapes, the dog is often so close that it looks like they are overlapping, but this is an abstract duck with no size, only a location, so unless the dog is exactly at the right spot, the duck escapes if it crosses the 1 unit border of the pond. Also, my strategy for the dog is to just go to the point on the edge of the pond that is closest to the duck and wait like a good boy, who’s a good boy, you’re a good boy. I’m not sure if that’s optimal but any movement from that spot seems like a mistake until the duck changes course or crosses the center point. This is one are where I think I am probably wrong.

The other area where I could be am probably wrong is in the way that I have the Smart duck scan the shore for a likely escape route. I just choose a range of shoreline and increment test locations in discrete chunks looking for a candidate. At the very least, a mathy solution would do that continuously, not discretely.

I wrote this app using MFC and Gdiplus so I’m not going to include the whole project here. I’ll just post the pond.h and pond.cpp files, which is where all the calculations and simulation are being done. Everything else is just presentation.

First off, I assume the pond has a radius of 1 unit with its center at 0,0, and that the duck moves at a rate of 1 unit/second. I store an x- and y-variable for the duck and start it at the center of the pond.

I could have used x- and y-variables for the dog, but since it always moves along the edge of the circular pond, I instead store its position as a single real value from -PI to PI. (It really should be 0 to TAU, because PI is wrong, no really, but whatevs.) I start the dog at 0, meaning the East edge of the pond.

The dog’s speed is expressed in terms of the duck’s speed (1/sec) * TAU/2 (okay, fine, PI.) So if the dog moves at PI * the duck’s speed, he can make a half orbit of the pond by the time the duck moves from the center to the edge. Using these basic units makes the code very simple since the time required for the duck to reach the edge is equal to 1 minus its distance from the center of the pond. The time required for the dog to reach any point of the pond’s edge is simply the angle to that point minus the dog’s position divided by the dog’s speed.

I also created a simple DPoint class with some convenient operators to make things a little easier. 

// Pond.h

#pragma once
#include <math.h>

#define PI 3.1415926535897932384626433832
#define TAU (PI * 2)

enum EDuck {
    DumbDuck,
    CowardlyDuck,
    CleverDuck,
    SmartDuck,
};

struct DPoint
{
    double  x;
    double  y;

    void   Set(double _x, double _y) {x = _x; y = _y;}
    double Magnitude() const
    {
        return sqrt(x * x + y * y);
    }
};
static DPoint operator -(const DPoint &l, const DPoint &r)
{
    DPoint pt;
    pt.x = l.x - r.x;
    pt.y = l.y - r.y;
    return pt;
}
static DPoint operator +(const DPoint &l, const DPoint &r)
{
    DPoint pt;
    pt.x = l.x + r.x;
    pt.y = l.y + r.y;
    return pt;
}
static DPoint operator *(const DPoint &l, const DPoint &r)
{
    DPoint pt;
    pt.x = l.x * r.x;
    pt.y = l.y * r.y;
    return pt;
}
static DPoint operator *(const DPoint &l, double scalar)
{
    DPoint pt;
    pt.x = l.x * scalar;
    pt.y = l.y * scalar;
    return pt;
}

class Pond;

class DuckStrategy
{
public:
    virtual DPoint  GetDestinationVector(Pond *pond) = 0;
};

class DumbDuckStrategy : public DuckStrategy
{
public:
    virtual DPoint  GetDestinationVector(Pond *)
    {
        DPoint pt;
        pt.x = -1.0;
        pt.y = 0;
        return pt;
    }
};

class CowardlyDuckStrategy : public DuckStrategy
{
public:
    virtual DPoint  GetDestinationVector(Pond *pond);
};

class CleverDuckStrategy : public CowardlyDuckStrategy
{
public:
    virtual DPoint  GetDestinationVector(Pond *pond);
};

class SmartDuckStrategy : public CowardlyDuckStrategy
{
public:
    virtual DPoint  GetDestinationVector(Pond *pond);
};

class Pond
{
    double  m_dogPos;   // in radians, 0-tau
    double  m_dogSpeed; // in radians/sec
    DPoint  m_duckPos;    
    DPoint  m_lastDuckVector;
    EDuck   m_behavior;
    bool    m_escaped;
    bool    m_beeline;
    double  m_elapsed;
    double  m_beelineTime;

    CowardlyDuckStrategy    m_cowardly;
    DuckStrategy    *m_strategy;

public:
    Pond() : m_dogPos(0)
    {
        SetStrategy(&m_cowardly);
        Restart(PI);
        m_escaped = true;   // so we have to hit restart to start
    }

    void    Restart(double dogSpeed)
    {
        m_beeline = false;
        m_elapsed = 0;
        m_beelineTime = 0;
        m_dogPos = 0;
        m_dogSpeed = dogSpeed;
        m_escaped = false;
        m_duckPos.Set(0, 0);
        m_lastDuckVector.Set(0, 0);
    }

    void    SetStrategy(DuckStrategy *strat)
    {
        m_strategy = strat;
    }

    void    SetBeeline() {m_beeline = true; m_beelineTime = m_elapsed;}
    bool    GetBeeline() const {return m_beeline;}
    double  GetElapsed() const {return m_elapsed;}
    DPoint  GetDuckPos() const {return m_duckPos;}
    double  GetDogPolar() const {return m_dogPos;}
    double  GetDogSpeed() const {return m_dogSpeed;}
    DPoint  GetLastDuckVector() const {return m_lastDuckVector;}
    DPoint  GetDogPos() const 
    {
        DPoint pos;
        pos.x = cos(m_dogPos);
        pos.y = sin(m_dogPos);
        return pos;
    }
    double  GetDogTimeToDest(double dest) const;
    double  GetDuckTimeToDest(double dest) const;

    void    Process(double sec)
    {
        if(m_escaped || m_strategy == NULL)
            return;


        // get duck vector
        DPoint duckVec = m_strategy->GetDestinationVector(this);
        m_lastDuckVector = duckVec;

        // get dog destination
        double mag = m_duckPos.Magnitude();
        double dogDest = PI;
        if(mag != 0.0)
        {
            double s = m_duckPos.y / mag;
            double c = m_duckPos.x / mag;
            dogDest = atan2(s, c);
        }

        // update duck position
        duckVec = duckVec * sec;
        m_duckPos = duckVec + m_duckPos;
        mag = m_duckPos.Magnitude();
        if(mag >= 1.0)
        {
            double overage = mag - 1.0;
            m_duckPos = m_duckPos * (1 / mag);
            sec -= overage;
            m_escaped = true;
        }
        m_elapsed += sec;
        if(mag < sec)
        {
            m_duckPos.Set(0, 0);
        }
        

        // update dog position
        double dogDelta = dogDest - m_dogPos;
        if(dogDelta < -TAU)
            dogDelta += TAU;
        else if(dogDelta > TAU)
            dogDelta -= TAU;

        double dogProgress = m_dogSpeed * sec;
        if(dogProgress > abs(dogDelta))
        {
            m_dogPos = dogDest;
        }
        else
        {
            m_dogPos += dogProgress;
            if(m_dogPos > PI)
                m_dogPos -= TAU;
            else if(m_dogPos < -PI)
                m_dogPos += TAU;
        }
    }
};

pond.h

// Pond.cpp

#include "stdafx.h"
#include "Pond.h"

DPoint CowardlyDuckStrategy::GetDestinationVector(Pond *pond)
{
    DPoint dogPos = pond->GetDogPos();
    DPoint duckPos = pond->GetDuckPos();

    DPoint pt = duckPos - dogPos;

    double mag = pt.Magnitude();
    if(mag == 0)
    {
        pt.x = -1.0;
        pt.y = 0;
    }
    else
    {
        pt.x /= mag;
        pt.y /= mag;
    }
    return pt;
}

DPoint CleverDuckStrategy::GetDestinationVector(Pond *pond)
{
    DPoint curVec = pond->GetDuckPos();
    double mag = curVec.Magnitude();

    if(mag == 0)
        return CowardlyDuckStrategy::GetDestinationVector(pond);

    double toGoal = 1.0 - mag;

    double dogPos = pond->GetDogPolar();
    double s = curVec.y / mag;
    double c = curVec.x / mag;
    double dogDest = atan2(s, c);

    double dogTime = pond->GetDogTimeToDest(dogDest);

    if(dogTime > toGoal)
    {
        pond->SetBeeline();
        curVec = curVec * (1 / mag);
        return curVec;
    }

    return CowardlyDuckStrategy::GetDestinationVector(pond);
}

DPoint SmartDuckStrategy::GetDestinationVector(Pond *pond)
{
    DPoint curVec = pond->GetDuckPos();
    double mag = curVec.Magnitude();

    if(pond->GetBeeline())
    {
        return pond->GetLastDuckVector();
    }

    if(mag == 0)
        return CowardlyDuckStrategy::GetDestinationVector(pond);

    double s = curVec.y / mag;
    double c = curVec.x / mag;
    double dest = atan2(s, c);

    double start = dest;
    double end = dest + PI / 8.0;
    double inc = (end - start) / 100;

    for(double test = start; test <= end; test += inc)
    {
        double duckTime = pond->GetDuckTimeToDest(test);
        double dogTime = pond->GetDogTimeToDest(test);
        if(duckTime < dogTime)
        {
            DPoint pt;
            pt.Set(cos(test), sin(test));

            pt = pt - curVec;

            pt = pt * (1 / pt.Magnitude());
            double one = pt.Magnitude();
            pond->SetBeeline();
            return pt;
        }
    }

    return CowardlyDuckStrategy::GetDestinationVector(pond);
}



double Pond::GetDuckTimeToDest(double dest) const
{
    DPoint pt;
    pt.x = cos(dest);
    pt.y = sin(dest);

    pt = pt - m_duckPos;
    return pt.Magnitude();
}

double Pond::GetDogTimeToDest(double dogDest) const
{
    double dogPos = GetDogPolar();

    double dogDelta = dogDest - dogPos;
    if(dogDelta < -PI)
        dogDelta += TAU;
    else if(dogDelta > PI)
        dogDelta -= TAU;

    double degrees = dogDelta * 360.0/PI;

    dogDelta = abs(dogDelta);
    ASSERT(dogDelta <= PI);

    double dogTime = dogDelta / GetDogSpeed();
    return dogTime;
}

pond.cpp

I chose to wrap most of the code for my pond in a single object that owns the duck and the dog. I could have chosen to make separate classes for the duck and dog and maybe have multiple ducks and multiple dogs each taking different approaches but then I thought, what if I didn’t instead.

As I said, I don’t think I got this one right, so I’ll be really interested to see what Ollie posts as the answer next Friday. Thanks for reading!