The Bunny Problem

Added 2021-12-04, Modified 2022-03-12

An approach using probability theory

The intended solution is some boring algorithm I'm not going to mention here, we're going to explore how to approach this from a probability theory standpoint (similar to Solving Wordle using information theory except my approach doesn't work)

We're going to keep a probability distribution over where the bunny could be, we start off with equal likelyhood for each hole

N = 100
P = [1/N] * N

Now we write the function to update our belief when we learn the bunny isn't in hole j

def guess(P, j):
  # after setting P[j] = 0 we have removed P[j] probability mass from
  # the distribution, so the new sum is 1-P[j]. Dividing everything by
  # 1-P[j] fixes this. (a minor detail, instead of 1 we use sum(P) for
  # numerical stability since sum(P) isn't exactly 1)
  s = sum(P) - P[j]
  return [
    0 if i == j else P[i] / s
    for i in range(N)

Now for the tricky bit, what happens to our beliefs when the bunny jumps to an adjacent hole?

It might seem hopelessly complicated, but you can figure it out by conditioning on three different cases

This is basically code, the one thing to be careful of is mutation- If we modified P[j] implace it would screw up the computation of P[j+1]. Hence we make a copy

def jumps(P):
  Pn = P.copy()
  Pn[0] = P[1]*0.5
  Pn[N-1] = P[N-2]*0.5

  for j in range(1, N-1):
    Pn[j] = 0.5*P[j-1] + 0.5*P[j+1]

  return Pn

We're pretty much done, we can wrap it up to test it

import random
def argmax(l):
  return max(range(len(l)), key=lambda i: l[i])

hole = random.randint(0, N-1)
t = 0
while True:
  t += 1
  j = argmax(P)
  if j == hole:
    print(f'Found the bunny in hole {j} after {t} tries!')
    print(f'[{t}] Guessed {j} with P[{j}] = {P[j]:.4f} (true is {hole})')
    # wrong hole, update our beliefs
    P = guess(P, j)
    P = jumps(P)

    # the bunny moves!
    if    hole == 0:   hole += 1
    elif  hole == N-1: hole -= 1
    else:              hole += random.choice([-1, 1])

Now we can test it

$ python3
[1] Guessed 0 with P[0] = 0.0100 (true is 93)
[2] Guessed 2 with P[2] = 0.0101 (true is 92)
[3] Guessed 4 with P[4] = 0.0102 (true is 93)
[42] Guessed 82 with P[82] = 0.0176 (true is 86)
[43] Guessed 84 with P[84] = 0.0179 (true is 87)
Found the bunny in hole 86 after 44 tries!

Cool, I wonder how well it does on average? If we factor the above loop into a function we can do something like

r = [play() for _ in range(10000)]

import statistics as stats
mean, median, stdev = stats.mean(r), stats.median(r), stats.stdev(r)
print(f'mean: {mean} median: {median} stddev {stdev}')

Running this I get

mean: 94.7903 median: 53.0 stddev 118.2896460809102

How does this compare to random guessing? replacing the j = argmax(P) line with j = random.randint(0, N-1) gives

mean: 99.5563 median: 70.0 stddev 97.57317399130908

Our algorithm barely beats random guessing :( something must have gone wrong, I'll figure it out later

The mean for random guessing being close to 100 makes sense since every time we guess randomly we have a p = 1/100 (independent) chance of success, which is exactly the Geometric Distribution meaning the expected value is 1/p = 100 lining up with our numerical result of around 99.6.

PS: Full code is downloadable here