## The Sultan’s Harem, Magic Parrot, and the Mathematician

Here’s an old problem. A mathematician is in a harem with 10 rooms. Each room has a different number of girls in it. The mathematician is blindfolded but has a magic parrot that knows probability. The mathematician needs to estimate the number of girls in each room. He can do so by using this algorithm:

```start in any room
loop max_iterations
parrot counts number girls in room
mathematician picks a random next room
parrot peeks into selected next room

if next room has more girls:
parrot says move to next room
else if next room has fewer girls:
move to next room with small p else stay

spend night in curr room and increment room count
end-loop

estimated dist of girls = room counts / iterations
```

In the demo program below, a harem is set up with 10 rooms where the first room has one girl, the second room has two girls, and so on. There are a total of 1 + 2 + 3 + . . + 10 = 55 girls in the harem. Therefore the true probability distribution to estimate is [1/55, 2/55, . . 10/55] = [0.0182, 0.0364, . . 0.1818].

The mathematician is estimating the probability distribution of girls in each room. There’s a ton of things going on here but the key idea is that it’s possible to estimate a probability distribution using comparisons of two states at a time. This is fascinating math related to Gibbs Sampling and Metropolis Sampling and Monte Carlo methods.

The clever part of the algorithm is when the parrot decides to tell the mathematician to move to a worse room (fewer girls). That happens with probability next_girls / curr_girls. For example, if the current room has 6 girls and the next room has 2 girls, the magic parrot will tell the mathematician to move to the worse room with p = 2 / 6 = 0.3333.

```# harem_parrot.py
# estimate a probability distribution

import numpy as np

def print_vec(v):
for i in range(len(v)):
print("%0.4f " % v[i], end="")
print("")

def main():
print("\nBegin sultan's harmem and magic parrot demo \n")
np.random.seed(0)

# set number  girls in each room
harem = np.array([1,2,3,4,5,6,7,8,9,10], dtype=np.float32)
harem /= np.sum(harem)  # a p-distribution

counts = np.zeros(10, dtype=np.int)  # nights spent each room
max_iter = 10000
curr_room = 5  # start here

for i in range(max_iter):
curr_girls = harem[curr_room]  # num girls curr room
next_room = np.random.randint(low=0, high=10)
next_girls = harem[next_room]

if next_girls > curr_girls:  # next room better; move
curr_room = next_room
else:   # next room is worse but move sometimes
t = next_girls / curr_girls   # threshhold
p = np.random.random() # [0.0, 1.0)
if p < t:  # accept and move to worse room
curr_room = next_room
else:  # reject and stay in same room
curr_room = curr_room  # for readability

counts[curr_room] += 1  # spend night in room

est_dist = counts / max_iter  # estimated distribution

print("True and estimateds distribution of girls: \n")
print_vec(harem)        # the true distribution
print_vec(est_dist)     # the estimated

print("\nEnd demo ")

if __name__ == "__main__":
main()
```

I don’t know much about harems, but I doubt that Hollywood depictions resemble reality. From left, “Son of Sinbad” (1955) with Vincent Price and the interesting Lili St. Cyr, “Don Juan DeMarco” (1994) with Johnny Depp, “Lost in a Harem” (1944) with Abbott and Costello, and “Harem Scarem” (1965) with Elvis Presley.