This post is a response to Zach Wissner-Gross’s Riddler puzzle on the shell game.

Riddler Express

Here is the puzzle description:

While walking down the street one day, you notice people gathering around a woman playing the shell game. With each game, she places a ball under one of three cups, and then swaps the positions of pairs of cups several times before asking passersby which cup they think the ball is now under.

You have it on good authority that she is playing fairly, performing all the moves in plain sight, albeit too fast for you to track precisely which cups she’s moving. However, you do have one additional key piece of information — every time she swaps cups, one of them has the ball. In other words, she never swaps the two empty cups.

When it’s your turn to guess, you note which cup she initially places the ball under. Then, as she begins to swap cups, you close your eyes and count the number of swaps. Once she is done, you open your eyes again. What is your best strategy for guessing which cup has the ball?

There are three cups, call them A, B and C. After any number of swaps, the probability that the ball is in one of those cups make a distribution. $P_{A,t} + P_{B,t} + P_{C,t}=1$.

At time $t$, we know that $P_{A,0}=1$. And, at every $t > 0$, the ball moves from its current position to one of the unoccupied cups.

$$P_{A,t} = (1/2) P_{B,t-1} + (1/2) P_{C, t-1} \quad P_{B,t} = (1/2) P_{A,t-1} + (1/2) P_{C, t-1} \quad P_{C,t} = (1/2) P_{A,t-1} + (1/2) P_{B, t-1} $$

This is a Markov Chain problem. The probabilty distribution at time $t$ is given by the following vector

$$p_t = \begin{pmatrix} 0 & 1/2 & 1/2 \newline 1/2 & 0 & 1/2 \newline 1/2 & 1/2 & 0 \end{pmatrix}^t \begin{pmatrix} 1 \newline 0 \newline 0\end{pmatrix}$$

I used Mathematica to exponentiate the matrix.

mat = {{0, 1, 1}, {1, 0, 1}, {0, 1, 1}}/2;
FullSimplify[MatrixPower[mat, t].Transpose[{{1, 0, 0}}], 
 Assumptions -> {n > 0}]

which gives me $$P_{A,t} = \frac{1}{6}+\frac{1}{3} \left(-\frac{1}{2}\right)^t \quad P_{B,t} = \frac{1}{6}-\frac{2}{3}\left(-\frac{1}{2}\right)^t \quad P_{C,t} = \frac{1}{6} + \frac{1}{3} \left(-\frac{1}{2}\right)^t$$

Note, that for $t$ odd, $P_{B,t} > P_{A,t}$ and $P_{B,t} >P_{C,t}$ and for $t$ even $P_{A,t}=P_{C,t} > P_{B,t}$

Therefore, the solution is to choose B if the number of swaps is odd, othewise A and C with equal probability.

Riddler Classic

Having solved the shell game in this week’s Express, you are now ready to play the quantum shell game. Instead of a ball, you are now trying to capture an electron. Now, you’re not sure precisely where it is, but you know it’s somewhere on a two-dimensional surface. What’s more, you know that the probability distribution is a 2D Gaussian (or “normal”) distribution. More precisely, the probability the electron is a distance r units in any direction from some central point is proportional to $exp(−r^2/2)$.

You have four cylindrical cups, each of which has a radius of 1 unit. In this game, you want to place the cups over the surface to maximize the probability that the electron is in one of the cups.

How should you place the cups, and what is the probability you will catch the electron?

Extra credit: What if you have different numbers of cups, like three, five, six or even more? How would you place the cups, and what would be your chances of catching the electron in one of them?

The probability mass of the electron is $\propto e^{-r^2/2}$. Let’s go ahead and normalize it such that it sums to one

$$ 4\pi \int_0^1 dr (k e^{-r^2/2}) = 1\implies k = \frac{1}{2\pi} \implies p(r) = \frac{1}{2\pi} e^{-r^2/2}$$

Let’s take a computational view of solving this problem. We are solving an optimization problem, with some contraints. The constraints being that the cups are physical objects that can’t overlap.

Consider a cup located at vector $\vec{v}$. The probabilty mass function it covers is given by $$\begin{align*} p_{\vec{v}} &= \frac{1}{2}\int_0^{2\pi} \int_0^1 e^{-r^2/2} \rho : \mathrm{d}\rho \mathrm{d}\theta \end{align*}$$

We have $r^2 = (\vec{v} - \vec{r}_{\rho, \theta})^2 = v^2 + r^2 - 2v r \cos{\theta}$.

This gives us $$p_{\vec{v}} = \frac{e^{-v^2/2}}{2\pi} \int_0^1 \mathrm{d}\rho : e^{-\rho^2} \mathrm{d}\rho \int_0^{2\pi} \mathrm{d}\theta e^{v\rho \cos{\theta}} \rho $$

We need to find 4 (or in general $n$) points such that the sum of $p_{\vec{v}_i}$ is maximized while making sure that the circles don’t overlap.

The maximization problem can be written as $$ \begin{align} \text{maximize } &\sum_{i=1}^{N} p_{\vec{v}_i} \newline \text{ s.t} \quad \quad &|\vec{v}_i-\vec{v}_j| >= 2 \quad \forall i,j : i \neq j \end{align} $$ The constraints ensure that the central locations of the cups are separated enough that their overlap is zero.

Here is the code that does this optimization.

import numpy as np
from scipy.integrate import dblquad
from scipy.optimize import minimize, NonlinearConstraint



def integral_under_cup(x, y):
    """
    Function to calculate the pmf under the cup situated at v_x, v_y
    """
    d2 = x**2 + y**2
    d = np.sqrt(d2)
    return np.exp(-d2/2)/(2*np.pi) * dblquad(lambda r, t: np.exp(-(r**2)/2)*r*np.exp(r*d*np.cos(t)), 
                                            0, 2*np.pi, 0, 1)[0]

def total_pmf_under_cups(x_s):
    """
    total pmf under all cups
    x_s: list of x and y coodiantes of N cups. 
    Even positions are x-coordinate.
    Odd positions are y-coordiante
    """
    acc = 0
    for i in range(len(x_s)//2):
        acc += integral_under_cup(x_s[2*i], x_s[2*i+1])
    return -acc

def no_overlap_constraint(x_s):
    """
    Returns a vector of separation between all pairs of cups.
    The optimize function should constraint it to be stricly greater than zero.
    """
    distance = []
    for i in range(len(x_s)//2):
        for j in range(i):
            x_a, y_a = x_s[2*i], x_s[2*i+1]
            x_b, y_b = x_s[2*j], x_s[2*j+1]
            distance.append((x_a - x_b)**2 + (y_a - y_b)**2 - 4)
            
    return distance

# define the constraint
constraints = NonlinearConstraint(no_overlap_constraint, 0, np.inf)


# Fix the problem size (number of cups)
N=4

# Now optimize!
optimal_solution = minimize(total_pmf_under_cups, x0=np.random.uniform(-1, 1, 2*N), constraints=[constraints])

# Here are your final positions
optimal_positions = optimal_solution.x

I did this for $N \in [4, 12]$. The results are shown below.

image