<jupyter_start><jupyter_text>Problem statementfrom : https://www.janestreet.com/puzzles/figurine-figuring/  Jane received 78 figurines as gifts this holiday season:  12 drummers drumming, 11 pipers piping, 10 lords a-leaping, etc., down to 1 partridge in a pear tree.   They are all mixed together in a big bag.  She agrees with her friend Alex that this seems like too many figurines for one person to have, so she decides to give some of her figurines to Alex.   Jane will uniformly randomly pull figurines out of the bag one at a time until she pulls out the partridge in a pear tree, and will give Alex all of the figurines she pulled out of the bag (except the partridge, that’s Jane’s favorite).    If n is the maximum number of any one type of ornament that Alex gets, what is the expected value of n, to seven significant figures?<jupyter_code>from math import factorial

import numpy as np
from tqdm.notebook import tqdm

mat_combis = [80 * [None] for _ in range(80)]


def n_combis(k, n):
    if mat_combis[k][n] is not None:
        return mat_combis[k][n]
    v = factorial(n) // (factorial(k) * factorial(n - k))
    mat_combis[k][n] = v
    return v


def compute_proba_n_drawn(n, n_drawn):
    p = 1
    n_tot = sum(range(1, n + 1))
    for i in range(n_drawn):
        p *= (n_tot - i - 1) / (n_tot - i)
    p *= 1 / (n_tot - n_drawn)
    return p


def rec(n, max_count, combi, n_drawn_above, n_drawn, n_below, max_filled):
    global s
    len_combi_max = n + 1 - max(2, max_count)
    n_drawn_below = n_drawn - n_drawn_above
    n_drawnable_max_future = sum(min(count, max_count) for count in range(2, n + 1 - len(combi)))
    if n_drawn_above > n_drawn or n_drawn_above + n_drawnable_max_future < n_drawn:
        return
    if len(combi) == len_combi_max:
        p = 1
        for count_tot, count_drawn_above in enumerate(combi, max(2, max_count)):
            p *= n_combis(count_drawn_above, count_tot)
        p *= n_combis(n_drawn_below, n_below)
        s += p
        return

    if len(combi) == len_combi_max - 1 and not max_filled:
        counts = [max_count]
    else:
        counts = range(max_count + 1)
    for count in counts:
        combi.append(count)
        rec(n, max_count, combi, n_drawn_above + count, n_drawn, n_below, max_filled or count == max_count)
        combi.pop()


def compute_exact_answer(n):
    global s
    esp_max_count = 0
    sum_proba_max_count = 0
    n_total = sum(range(2, n + 1))
    for max_count in tqdm(range(0, n + 1)):
        proba_max_count = 0
        n_below = sum(c for c in range(2, n + 1) if c < max_count)
        for n_drawn in tqdm(range(max_count, sum(range(2, n + 1)) + 1)):
            s = 0
            rec(n, max_count, [], 0, n_drawn, n_below, max_filled=False)
            n_permuts_drawn_max_count_ok = s * factorial(n_drawn)
            n_permuts_drawn_total = factorial(n_total) / factorial(n_total - n_drawn)
            proba_max_count_given_n_drawn = n_permuts_drawn_max_count_ok / n_permuts_drawn_total
            proba_max_count += proba_max_count_given_n_drawn * compute_proba_n_drawn(n, n_drawn)
        sum_proba_max_count += proba_max_count
        esp_max_count += proba_max_count * max_count
    assert abs(sum_proba_max_count - 1) < 1e-10
    return esp_max_count


def approximate_answer_with_monte_carlo(n, n_iter=2_000_000):
    l = []
    for i in range(1, n + 1):
        l.extend(i * [i])
    list_n = []
    for _ in tqdm(range(n_iter)):
        np.random.shuffle(l)
        index_1 = l.index(1)
        list_n.append(index_1 and np.bincount(l[:index_1]).max())
    return np.mean(list_n)<jupyter_output><empty_output><jupyter_text>Checking code against Monte Carlo simulations<jupyter_code>%%time
for N in range(2,11):
    print("N =", N)
    exact = compute_exact_answer(N)
    approx = approximate_answer_with_monte_carlo(N)
    error = abs(exact-approx)
    print(N, exact, approx, error)
    assert error < 0.005<jupyter_output>N = 2<jupyter_text>Computing the solution to the problem<jupyter_code>%%time
answer = compute_exact_answer(12)
answer<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.093994140625, 'fraction_numerical': 0.016357421875, 'mean_word_length': 3.2150205761316872, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 13, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '549571', 'license': 'MIT', 'n_tokens_mistral': 1582, 'n_tokens_neox': 1467, 'n_words': 415, 'path': 'JaneStreet_puzzle_jan_2021_Figurine_Figuring.ipynb', 'repo_name': 'ismael-elatifi/jane-street-puzzles', 'size': '4096'}