1 January 2022 · Problems Maths

Robot Archery

Original problem

Happy New Year!

Jane Street had a doozy of a December problem, which is linked above. It had some similarity to the August problem, involving robots, probability and differential equations.

In this month's archery contest, four identical robots take turn to strike a (large finite) target at random. One key insight here is that the size and distribution of the shots don't matter as we're only concerned about the orderings of those shots (ie. it doesn't matter how far the shots are spread from each other; all that matters is if a certain shot is closer or further to the bullseye than another). This means that without loss of generality, we can assume that shots (more specifically, distance from the bullseye) are spread uniformly along $[0,1]$, where 0 is the bullseye and 1 (in some units) is the furthest point from the bullseye.

Define $$f_{i,j}(x)$$as the probability that if there are $i$ total robots in the queue, the $j$th robot in the queue eventually wins if the current marker (ie. closest shot so far) is at $x$.

So, we want to find $f_{4,4}(1)$. Let's first scale down and consider a competition with only two robots.

Two robots

Consider that both robots are still in queue and the marker is at $x$. What is the probability that Robot #1 loses (ie. $1-f_{2,1}(x)$)? Two things can happen that will result in Robot #1 losing:

This means we have $$\begin{align*}1-f_{2,1}(x) &= (1-x)+\int_0^{x}f_{2,1}(a)da\\\therefore f_{2,1}(x)&=x - \int_0^{x}f_{2,1}(a)da\end{align*}$$

Differentiating gives $$ f_{2,1}'(x)=1 - f_{2,1}(x)$$

With the obvious boundary condition of $f_{2,1}(0)=0$, the solution is $$\boxed{f_{2,1}(x)=1-e^{-x}}$$and so $$\boxed{f_{2,2}(x)=e^{-x}}$$Nice. Let's add another competitor.

Three robots

The logic for three robots is the same, but the maths becomes more tedious. Once again, consider the probability that Robot #1 loses if the marker is at $x$ currently. Here are the cases that eventually lead to Robot #1's loss:

The equation that accounts for all four cases (respectively) above is:

$$\begin{align*}1-f_{3,1}(x) &= (1-x)\\&\;\;\;+\int_0^{x}f_{3,1}(a)da\\&\;\;\;+\int_0^x \int_0^a f_{3,1}(b)dbda\\&\;\;\;+\int_0^x(1-a) f_{2,1}(a)da\end{align*}$$


$$\begin{align*}0 &= -x + f_{3,1}(x)\\&\;\;\;+\int_0^{x}f_{3,1}(a)da\\&\;\;\;+\int_0^x \int_0^a f_{3,1}(b)dbda\\&\;\;\;+\int_0^x(1-a) f_{2,1}(a)da\end{align*}$$

Differentiating both sides: $$\begin{align*}0 &= -1+f_{3,1}'(x) + f_{3,1}(x) + \int_0^x f_{3,1}(b)db + (1-x)f_{2,1}(x)\end{align*}$$

Before progressing, note that substituting $x=0$ to the above equation yields $f_{3,1}'(0)=1$ as an initial value. Now, differentiate again.

$$\begin{align*}0 &= f_{3,1}''(x) + f_{3,1}'(x) + f_{3,1}(x) + \frac{d}{dx}\Big[ (1-x)f_{2,1}(x)\Big]\end{align*}$$

Knowing $f_{2,1}(x) = (1-e^{-x})$, the above becomes

$$f_{3,1}''(x) + f_{3,1}'(x) + f_{3,1}(x) = xe^{-x} - 2e^{-x}+1$$

Applying Laplace transforms to both sides, and using the initial values $f_{3,1}(0)=0$ and $f_{3,1}'(0)=1$ gives:

$$\begin{align*}(s^2+s+1)\mathcal{L}(f_{3,1})(s) - 1 &= \frac{1}{(s+1)^2} - \frac{2}{s+1} + \frac{1}{s}\\\therefore \mathcal{L}(f_{3,1})(s) &=\frac{1}{s^2+s+1}\Big[\frac{1}{(s+1)^2} - \frac{2}{s+1} + \frac{1}{s}+1\Big]\end{align*}$$

The inverse transform is painful to do by hand (a lot of factorising!), but WolframAlpha comes to the rescue and reports the inverse transform as

$$\boxed{f_{3,1}(x)=1+e^{-x}(x-1) -e^{-x/2}\Big[\frac{2}{\sqrt{3}}\sin\Big(\frac{\sqrt{3}}{2}x\Big)\Big]}$$

Next, let's find $f_{3,2}(x)$. Given some $x$, Robot #2 wins if either the following occurs:

The equation hence is $$f_{3,2}(x)=(1-x)f_{2,1}(x) + \int_0^x f_{3,1}(a)da$$

Knowing both $f_{2,1}(x)$ and $f_{3,1}(x)$, the above evaluates to

$$\boxed{f_{3,2}(x)=-e^{-x} + e^{-x/2}\Big[\frac{1}{\sqrt{3}}\sin\Big(\frac{\sqrt{3}}{2}x\Big)+\cos\Big(\frac{\sqrt{3}}{2}x\Big)\Big]}$$

$f_{3,3}(x)$ is simply the complement of $f_{3,1}(x)$ and $f_{3,2}(x)$ so I won't bother writing it down here (plus, we won't need it).

It's time to add the fourth robot.

Four robots

Once again, the process is the same. We consider the number of ways that Robot #1 can lose:

These cases respectively contribute to the following equation:

$$\begin{align*}1-f_{4,1}(x)&=(1-x) + \int_0^xf_{4,1}(a)da\\&\;\;\;+\int_0^x\int_0^af_{4,1}(b)dbda\\&\;\;\;+\int_0^x\int_0^a\int_0^bf_{4,1}(c)dcdbda\\&\;\;\;+\int_0^x\int_0^a(1-b)f_{3,1}(b)dbda\\&\;\;\;+\int_0^x(1-a)f_{3,1}(a)da\\&\;\;\;+\int_0^x(1-a)\int_0^af_{3,1}(c)dcda\\&\;\;\;+\int_0^x(1-a)^2f_{2,1}(a)da\end{align*}$$

Differentiating and rearranging eventually results in this third-order ODE:


Note: at every differentiation step, we must evaluate the expression at $x=0$ to obtain the initial values:


To solve the ODE with the initial values, we again use Laplace transforms, WolframAlpha and a lot of paper later. This gives: $$\boxed{\begin{align*}f_{4,1}(x)&=1+e^{-x}\Big[\frac{-3}{4}x^2+2x-\frac{5}{4}\Big]\\&\;\;\;+e^{-x/2}\Big[\frac{3x-1}{\sqrt{3}}\sin\Big(\frac{\sqrt{3}}{2}x\Big)-(x+1)\cos\Big(\frac{\sqrt{3}}{2}x\Big)\Big]\\&\;\;\;+\frac{5}{4}\Big[\cos(x)-\sin(x)\Big]\end{align*}}$$

For Robot #2, the win-probability (which I won't explicitly write out here) can be determined by: $$f_{4,2}(x)=(1-x)f_{3,1}(x) + \int_0^x f_{4,1}(a)da$$

Once that's resolved, Robot #3 has $$f_{4,3}(x)=(1-x)f_{3,2}(x) + \int_0^x f_{4,2}(a)da$$

Finally, we arrive at $f_{4,4}(x)$, which is the complement of the other win probabilities. $$\boxed{\begin{align*}f_{4,4}(x)&=e^{-x}\Big[\frac{3}{4}x^2-\frac{7}{2}x+\frac{13}{4}\Big]\\&\;\;\;+e^{-x/2}\Big[\frac{5}{\sqrt{3}}\sin\Big(\frac{\sqrt{3}}{2}x\Big)+(2x-1)\cos\Big(\frac{\sqrt{3}}{2}x\Big)\Big]\\&\;\;\;-\frac{5}{4}\Big[\sin(x)+\cos(x)\Big]\end{align*}}$$

And so, the final answer to this problem is:

$$\color{red}{\boxed{\begin{align*}f_{4,1}(1)&=\frac{1}{2e}+\frac{1}{\sqrt{e}}\Big[\frac{5}{\sqrt{3}}\sin\Big(\frac{\sqrt{3}}{2}\Big)+\cos\Big(\frac{\sqrt{3}}{2}\Big)\Big]-\frac{5}{4}\Big[\sin(1)+\cos(1)\Big]\\&\approx0.1834376509\end{align*}}}$$to the required 10 decimal places.

The following code simulates this game $10^8$ times, which supports the above analytic result.

import random

N = 100000000
x = 0
for trial in range(N):
  players = ["A", "B", "C", "D"]
  marker = 1

  while len(players) > 1:
    throw = random.random()
    if throw < marker:
      marker = throw
      players = players[1:] + [players[0]]
      players = players[1:]
  if players[0] == "D":
    x += 1

# outputs 0.18348236