Return to Question

Edit:

Thanks a lot for your help @Peilonrayz and @Marc! I have modified my code according to your suggestions, and now it looks something like this:

import random
import statistics

def simulation(days_in_year):
    days = [0] * days_in_year
    num_emperors = 0
    holidays = 0
    #increments (as suggested) is slightly faster than while not all(days)
    while holidays < days_in_year:
        #also, slightly faster than np.random.randint. But why?
        bday = random.randint(0, days_in_year - 1)
        num_emperors += 1

        if days[bday] == 0:
            days[bday] = 1
            holidays += 1
        else:
            continue

        yesterday = (bday - 1) % days_in_year
        day_before_yesterday = (bday - 2) % days_in_year

        if days[day_before_yesterday] == 1 and days[yesterday] == 0:
            days[yesterday] = 1
            holidays += 1

        tomorrow = (bday + 1) % days_in_year
        day_after_tomorrow = (bday + 2) % days_in_year

        if days[day_after_tomorrow] == 1 and days[tomorrow] == 0:
            days[tomorrow] = 1
            holidays += 1
    return num_emperors


def monte_carlo(iters, days):
    for _ in range(iters):
        yield simulation(days)

The time required has become as followed:

avg_n_emperor = statistics.mean(monte_carlo(iters = 1000, days = 360)) #about 1.7 seconds

avg_n_emperor = statistics.mean(monte_carlo(iters = 100000, days = 360)) #about 170 seconds

avg_n_emperor = statistics.mean(monte_carlo(iters = 100000, days = 730)) #about 400 seconds

The correlation between time and iters & days is already becoming linear. I will try to modify a bit more and see where it goes :)

Edit:

Thanks a lot for your help @Peilonrayz and @Marc! I have modified my code according to your suggestions, and now it looks something like this:

import random
import statistics

def simulation(days_in_year):
    days = [0] * days_in_year
    num_emperors = 0
    holidays = 0
    #increments (as suggested) is slightly faster than while not all(days)
    while holidays < days_in_year:
        #also, slightly faster than np.random.randint. But why?
        bday = random.randint(0, days_in_year - 1)
        num_emperors += 1

        if days[bday] == 0:
            days[bday] = 1
            holidays += 1
        else:
            continue

        yesterday = (bday - 1) % days_in_year
        day_before_yesterday = (bday - 2) % days_in_year

        if days[day_before_yesterday] == 1 and days[yesterday] == 0:
            days[yesterday] = 1
            holidays += 1

        tomorrow = (bday + 1) % days_in_year
        day_after_tomorrow = (bday + 2) % days_in_year

        if days[day_after_tomorrow] == 1 and days[tomorrow] == 0:
            days[tomorrow] = 1
            holidays += 1
    return num_emperors


def monte_carlo(iters, days):
    for _ in range(iters):
        yield simulation(days)

The time required has become as followed:

avg_n_emperor = statistics.mean(monte_carlo(iters = 1000, days = 360)) #about 1.7 seconds

avg_n_emperor = statistics.mean(monte_carlo(iters = 100000, days = 360)) #about 170 seconds

avg_n_emperor = statistics.mean(monte_carlo(iters = 100000, days = 730)) #about 400 seconds

The correlation between time and iters & days is already becoming linear. I will try to modify a bit more and see where it goes :)

On planet J, a year lasts for D days. Holidays are defined by the two following rules.

1. At the beginning of the reign of the current Emperor, his birthday is declared a holiday from that year onwards.
2. If both the day before and after a day d are holidays, then d also becomes a holiday.

Initially there are no holidays. Let E(D) be the expected number of Emperors to reign before all the days of the year are holidays, assuming that their birthdays are independent and uniformly distributed throughout the D days of the year.

You are given E(2) = 1, E(5) = 31/6, E(365) is about 1174.3501.

Find E(1000) . Give your answer rounded to 4 digits after the decimal point.