empty-birthdays-in-december.py

#!/usr/bin/env python3

import decimal
import math

from decimal import Decimal
from fractions import Fraction

# In response to https://www.reddit.com/r/mildlyinteresting/comments/1hcb9ce/not_a_single_person_at_my_2000_student_high/
# and, in particular, some sloppy math in https://www.reddit.com/r/mildlyinteresting/comments/1hcb9ce/comment/m1mzfgh/?utm_source=share&utm_medium=web3x&utm_name=web3xcss&utm_term=1&utm_content=share_button
#
# Note that non-uniformity likely plays a very large role in the real distribution.
# This is based on a data-free prior.

# independent events: 0.120689928872311369108421247197 0.879310071127688630891578752803
# no leap years:      0.120795975400973612354880396771 0.879204024599026387645119603229
# with leap yars:     0.121222615424473694930263066452 0.878777384575526305069736933548

# Note that the leap-year accounting has a larger impact on the final probabilities
# than removing the independence relaxation.

def main():
    # Most naively assume that each day is independent of the others (reasonable
    # assumption when the individual probabilities are low).
    p = compute_independent_prob(Fraction(1, 365), 2000)
    print("independent events:", render_frac(p, 30), render_frac(1 - p, 30))
    # Naively assume that each day in the year has 1/365 uniform probability.
    p = compute_prob(Fraction(1, 365), 2000)
    print("no leap years:     ", render_frac(p, 30), render_frac(1 - p, 30))
    # Slightly less-naively account for leap years. In a leap year, the naive
    # probability of a birth on any given day in December is 1/366. In a
    # non-leap year, the probability is 1/365. Leap years occur _approximately_
    # every 1 of 4 years or, more precisely, every 97 of 400 years. Therefore,
    # the total probability of any given day in December is
    # 97/400 * 1/366 + 303/400 * 1/365. Note that this will actually give _less_
    # accurate results for present day, since the year 2000 _was_ a leap year
    # and we can safely assume there arent many people around from 1900 or
    # earlier. (On top of this, the school likely only includes birthdays from
    # a very small span of years.) For that reason, we stick to 1/4.
    p = compute_prob(
        Fraction(1, 4) * Fraction(1, 366) + Fraction(3, 4) * Fraction(1, 365),
        2000)
    print("with leap yars:    ", render_frac(p, 30), render_frac(1 - p, 30))


def compute_independent_prob(daily_prob: Fraction, students: int) -> Fraction:
    days = 31
    p_day_empty = (1 - daily_prob)**students
    p_approx_no_empty_days = (1 - p_day_empty)**31
    return 1 - p_approx_no_empty_days


def compute_prob(daily_prob: Fraction, students: int) -> Fraction:
    days = 31
    total_p = Fraction(0, 1)
    for i in range(1, days + 1):
        c = math.comb(days, i)
        # NOTE: The probability of a single birthday happening on any given day
        # is mutually exclusive with other days.
        p = c * (1 - i * daily_prob)**students
        if i % 2 == 0:
            total_p -= p
        else:
            total_p += p
    return total_p


def render_frac(frac: Fraction, decimal_places: int) -> str:
    scale = 10**decimal_places
    whole = whole_digits(frac)
    with decimal.localcontext() as ctx:
        # NOTE: We give ourselves 4 extra digits of wiggle room for precision.
        # The actual format string is computed below.
        ctx.prec = whole + decimal_places + 4
        rounded = Decimal(round(frac * scale)) / Decimal(scale)
        fmt = f"1e-{decimal_places}"
        return str(rounded.quantize(Decimal(fmt)))


def whole_digits(frac: Fraction) -> int:
    if frac < 0:
        raise Exception("fractions must be non-negative")
    whole_part = math.floor(frac)
    return ilog10(whole_part)


def ilog10(n: int) -> int:
    m = 1
    digits = 0
    while n >= m:
        digits += 1
        m *= 10
    return digits


if __name__ == "__main__":
    main()

~LICENSE

MIT License

Copyright (c) 2024 Benjamin Sidhom

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.