# The Bunny Problem

Created 2021-12-04, Edited 2022-03-12
An approach using probability theory

• There are 100 holes in a line, and there is a bunny hiding in one of the holes
• You can only look in one hole at a time, and every time you look the bunny randomly jumps to an adjacent hole

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

• If j = 0 then the only way for the bunny to be in hole 0 is for it to be in hole 1 then move to hole 0, in other words P = P*0.5
• If j = N-1 then P[N-1] = P[N-2]*0.5 by the same logic
• Otherwise P[j] = P[j-1]*0.5 + P[j+1]*0.5 since we could get to hole j from the left or the right.

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 = P*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!')
break
else:
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 bunny.py
 Guessed 0 with P = 0.0100 (true is 93)
 Guessed 2 with P = 0.0101 (true is 92)
 Guessed 4 with P = 0.0102 (true is 93)
...
 Guessed 82 with P = 0.0176 (true is 86)
 Guessed 84 with P = 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.