Birthday paradox

import matplotlib.pyplot as plt
import numpy as np
import random
import subprocess

def ratio_bdays(n):
    mult = 1
    for i in range(365, 365 - n , -1):
        mult = mult * i / 365
    return (1 - mult) * 100

def create_matrix_ratio(num_of_sims, num_days_random):
    sim_with_coinc = 0

    #Creates a 12x31 blue matrix
    list_calendar = np.zeros((12,31))
    for i in range(0,12):
        for j in range(0,31):
            list_calendar[i][j] = 0

    #Creates the non-existing days at the end of the month
    ##Months with 30 days
    for i in [3, 5, 8, 10]:
        list_calendar[i][30] = 100
    #February
    for i in [28, 29, 30]:
        list_calendar[1][i] = 100

    #Creates random bdays and draw them and the coincidences
    sum_days_random = 0
    coincidences = 0

    while sum_days_random < num_days_random:
        i = random.randint(0,11)
        j = random.randint(0,30)

        if list_calendar[i][j] == 30:  #If there is already a birthday (orange) -> green
            list_calendar[i][j] = 45
            coincidences = coincidences + 1
            sum_days_random = sum_days_random + 1

        elif list_calendar[i][j] != 100: #If day exists and not a bday already
            list_calendar[i][j] = 30
            sum_days_random = sum_days_random + 1

    if coincidences != 0:
        sim_with_coinc = sim_with_coinc + 1

    ax = plt.gca()

    for i in range(0,11):
        ax.axhline(y= i + 0.45, linewidth=0.5, color='k')
    for i in range(0,30):
        ax.axvline(x= i + 0.45, linewidth=0.5, color='k')

    ax.set_title("Random birthdays (orange): " + str(num_days_random) + "\n"
                "Coincidences / Shared birthdays (green): " +
                str(coincidences) + "\n" +
                "Simulation number " + str(num_of_sims) + "\n")

    return list_calendar, num_of_sims, sim_with_coinc, coincidences


def creates_graph(simus_user, num_days_random):
    list_simu_0 = []
    list_simu_1 = []

    list_ratio_0 = []
    list_ratio_sum_0 = []

    list_ratio_1 = []
    list_ratio_sum_1 = []

    list_num_coinc_0 = []
    list_num_coinc_1 = []

    for num_of_sims in range(1, simus_user):
        #Creates and saves the graph
        plt.figure(figsize=(12,8))

        ###### First matrix [0,0]###
        ax1 = plt.subplot(221)
        list_calendar, simu_with_coinc_0, ratio_0, num_coinc_0 = create_matrix_ratio(num_of_sims, num_days_random[0])
        list_simu_0.append( simu_with_coinc_0 )
        list_ratio_0.append( ratio_0 )
        list_ratio_sum_0.append( 100 * sum(list_ratio_0) / num_of_sims )
        list_num_coinc_0.append(num_coinc_0)

        #Turns the list into an array and colours it
        array_calendar = np.asarray(list_calendar)
        array_calendar = array_calendar.reshape((12, 31))
        ax1.matshow(array_calendar, cmap="tab20c", interpolation='none')

        col_labels = range(1,32)
        ax1.set_xticks(range(31))
        ax1.set_xticklabels(col_labels, fontsize=7)
        row_labels = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dic" ]
        ax1.set_yticks(range(12))
        ax1.set_yticklabels(row_labels, fontsize=8)

        ###### First ratio [1,0]###
        ax1 = plt.subplot(223)
        ratio_obj = round( ratio_bdays( num_days_random[0]) , 1)
        ax1.axhline(y= ratio_obj, linewidth=0.5, color='g')
        ax1.set_xlim( [0, simus_user] )
        plt.title("Ratio of simulations with shared birthdays (%)" + "\n"
                        + "Expected ratio: " + str(ratio_obj) + "%" + "\n"
                        + "Average of shared birthdays: " +
                        str( round(np.mean(list_num_coinc_0), 2) ) + "\n")
        ax1.plot(list_simu_0, list_ratio_sum_0, "-ro", linewidth=1, markersize=2)


        ###### Second matrix [0,1]###
        ax1 = plt.subplot(222)

        list_calendar, simu_with_coinc_1, ratio_1, num_coinc_1 = create_matrix_ratio(num_of_sims, num_days_random[1])
        list_simu_1.append( simu_with_coinc_1 )
        list_ratio_1.append( ratio_1 )
        list_ratio_sum_1.append( 100 * sum(list_ratio_1) / num_of_sims )
        list_num_coinc_1.append(num_coinc_1)

        #Turns the list into an array and colours it
        array_calendar = np.asarray(list_calendar)
        array_calendar = array_calendar.reshape((12, 31))
        ax1.matshow(array_calendar, cmap="tab20c", interpolation='none')

        col_labels = range(1,32)
        ax1.set_xticks(range(31))
        ax1.set_xticklabels(col_labels, fontsize=7)
        row_labels = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dic" ]
        ax1.set_yticks(range(12))
        ax1.set_yticklabels(row_labels, fontsize=8)


        ###### Second ratio [1,1]###
        plt.subplot(224)
        ratio_obj = round( ratio_bdays( num_days_random[1]) , 1)
        plt.axhline(y= ratio_obj, linewidth=0.5, color='g')
        plt.xlim( [0, simus_user] )
        plt.title("Ratio of simulations with shared birthdays (%)" + "\n"
                        + "Expected ratio: " + str(ratio_obj) + "%" + "\n"
                        + "Average of shared birthdays: " +
                        str( round(np.mean(list_num_coinc_1), 2) ) + "\n")
        plt.plot(list_simu_1, list_ratio_sum_1, "-ro", linewidth=1, markersize=2)

        if len(str(num_of_sims)) == 1:
            zeros = "00"
        elif len(str(num_of_sims)) == 2:
            zeros = "0"
        else:
            zeros = ""

        plt.savefig(zeros + str(num_of_sims) + ".png", bbox_inches='tight')
        plt.close()
        plt.close('all')


number_of_sims_to_run = 150
random_days = [23, 41]
creates_graph(number_of_sims_to_run, random_days)

delay = 6
subprocess.run(["convert", "-delay", str(delay), "-loop", "0", "*.png", "ZZZ.gif"])