Python Puzzles
You die and arrive at the pearly gates of heaven. Unfortunately, the pearly gates are a set of three random doors. One door leads directly to heaven, one door leads to a one-day stay in purgatory, and one door leads to a two-day stay in purgatory. Once you complete your stay in purgatory, you are back again to the pearly gates, where the doors are again randomized.
What is the average number of days spent in purgatory?
On average, a person will spend three days in purgatory.
Upon running the simulation multiple times, I observed that the maximum number of days it takes to find heaven is typically between 50 and 60 days. It’s interesting to note that the probability of spending one day in purgatory is lower than that of spending two days (11% vs. 14%). This makes sense since one may choose the one-day door multiple times during the decision-making process, leading to two, three, or more days spent in purgatory.
By overlaying a geometric distribution on the resulting counts by day, we can confirm that the data follows a geometric distribution. This distribution models the number of trials it takes to achieve the first success in a sequence of independent trials with a fixed probability of success. The overlaid distribution provides a good fit to the observed data, further supporting the use of a geometric distribution to model this scenario.
Here are the counts per day.
0 day(s): 333522
1 day(s): 110539
2 day(s): 148830
3 day(s): 86842
4 day(s): 77695
5 day(s): 54902
6 day(s): 44506
7 day(s): 32961
8 day(s): 25715
9 day(s): 19407
10 day(s): 14950
11 day(s): 11549
12 day(s): 8891
13 day(s): 6969
14 day(s): 5340
15 day(s): 3971
16 day(s): 3116
17 day(s): 2431
18 day(s): 1873
19 day(s): 1395
20 day(s): 1051
21 day(s): 837
22 day(s): 618
23 day(s): 486
24 day(s): 373
25 day(s): 293
26 day(s): 213
27 day(s): 177
28 day(s): 131
29 day(s): 121
30 day(s): 81
31 day(s): 49
32 day(s): 33
33 day(s): 30
34 day(s): 26
35 day(s): 12
36 day(s): 11
37 day(s): 13
38 day(s): 10
39 day(s): 2
40 day(s): 9
41 day(s): 4
42 day(s): 2
43 day(s): 4
44 day(s): 2
45 day(s): 1
46 day(s): 2
47 day(s): 0
48 day(s): 1
49 day(s): 0
50 day(s): 1
51 day(s): 0
52 day(s): 0
53 day(s): 1
54 day(s): 1
55 day(s): 0
56 day(s): 0
57 day(s): 0
58 day(s): 1To keep it organized, the Python code is divided into separate sections for simulation, print output, and graphing. I hard-coded the output from the simulation into the graphing section for convenience, but it should be avoided whenever possible for better flexibility and modifiability.
import random
total_time = []
result_counts = {}
# run the simulation 1 million times
for _ in range(1_000_000):
time = 0
door = random.randint(1, 3)
# keep opening doors until door 1 is opened (heaven is found)
while door != 1:
# if door 2 is opened, add 1 day to the total time
if door == 2:
time += 1
# if door 3 is opened, add 2 days to the total time
elif door == 3:
time += 2
door = random.randint(1, 3)
# add the total time to the result list and increment the corresponding dictionary key
total_time.append(time)
if time in result_counts:
result_counts[time] += 1
else:
result_counts[time] = 1
# get the minimum and maximum value in the result list
min_time = min(total_time)
max_time = max(total_time)
# update the result_counts dictionary to include all values between min_time and max_time
for time in range(min_time, max_time+1):
if time not in result_counts:
result_counts[time] = 0
# print the average time to find the prize
print("Average time to find heaven:", sum(total_time) / len(total_time))
# print the count of each number of days
print("Counts:")
for time, count in sorted(result_counts.items()):
print(f"{time} day(s): {count}")
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import geom
data = {0: 333522, 1: 110539, 2: 148830, 3: 86842, 4: 77695, 5: 54902, 6: 44506, 7: 32961, 8: 25715, 9: 19407, 10: 14950, 11: 11549, 12: 8891, 13: 6969, 14: 5340, 15: 3971, 16: 3116, 17: 2431, 18: 1873, 19: 1395, 20: 1051, 21: 837, 22: 618, 23: 486, 24: 373, 25: 293, 26: 213, 27: 177, 28: 131, 29: 121, 30: 81, 31: 49, 32: 33, 33: 30, 34: 26, 35: 12, 36: 11, 37: 13, 38: 10, 39: 2, 40: 9, 41: 4, 42: 2, 43: 4, 44: 2, 45: 1, 46: 2, 47: 0, 48: 1, 49: 0, 50: 1, 51: 0, 52: 0, 53: 1, 54: 1, 55: 0, 56: 0, 57: 0, 58: 1}
# Convert the data to a numpy array for convenience
x, y = np.array(list(data.keys())), np.array(list(data.values()))
# Normalize the counts to form a probability mass function (PMF)
pmf = y / np.sum(y)
# Plot the PMF as a bar chart
fig, ax = plt.subplots()
ax.bar(x, pmf, label="Discrete Probability Distribution", alpha=0.5, width=0.5)
# Plot the overlaid geometric distribution with p = 1/3
geom_dist = geom(1/3)
x_geom = np.arange(geom_dist.ppf(0.01), geom_dist.ppf(0.99))
ax.plot(x_geom, geom_dist.pmf(x_geom), label="Geometric Distribution (p = 1/3)", color="red")
# Set the plot labels and legend
ax.set_xlabel("Number of days to find Heaven")
ax.set_ylabel("Probability")
ax.set_title("Probability Distribution")
ax.legend()
plt.savefig('pearly-gates-probability-distribution-graph.png', dpi=300, bbox_inches='tight')
plt.show()