Opening hand analysis

deck_draw_mc.py

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

opening_hand_size = 5
deck_sizes = range(40, 62, 2)
n_trials = 10000
rates = []
gt_rates = []
std = []
n_card_copies = 3
# You know, I thought using strings here would be slower.
# But I tried replacing them with ints, and it didn't get
# much faster. So the joke gets to live on.
desired_card_one = "Left Hand of Exodia"
desired_card_two = "Right Hand of Exodia"
trash_card = "Kuriboh"

for deck_size in deck_sizes:

    # Theoretical expectation:
    from math import factorial
    def choose(n, r):
        # https://stackoverflow.com/questions/3025162/statistics-combinations-in-python/3027128
        return float(factorial(n)) // factorial(r) // factorial(n-r)
    def hypergeometric_distrib(num_c1, num_c2, hand_size, deck_size):
        numerator = (choose(n_card_copies, num_c1) *
                     choose(n_card_copies, num_c2) *
                     choose(deck_size - n_card_copies - n_card_copies, hand_size - num_c1 - num_c2))
        return numerator / choose(deck_size, hand_size)
    total_expected = 0.
    for num_c1 in range(1, min(n_card_copies, opening_hand_size)):
        for num_c2 in range(1, min(n_card_copies, opening_hand_size-num_c1)):
            total_expected += hypergeometric_distrib(num_c1, num_c2, opening_hand_size, deck_size)
    gt_rates.append(total_expected)

    deck = [desired_card_one]*n_card_copies \
         + [desired_card_two]*n_card_copies \
         + [trash_card]*(deck_size-n_card_copies*2)
    outcomes = []
    for trial_k in range(n_trials):
        hand = np.random.choice(deck, size=opening_hand_size, replace=False)
        #print("Hand: ", ", ".join(hand))
        if desired_card_one in hand and desired_card_two in hand:
            outcomes.append(True)
        else:
            outcomes.append(False)

    rate = float(sum(outcomes))/float(n_trials)
    rates.append(rate)
    #print("Of %d trials, %d (= %f%%) opening hands contained both %s and %s." % (
    #    n_trials, sum(outcomes), rate, desired_card_one, desired_card_two))

# Normal approx https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval
z = 1.96
errorbars = z*np.sqrt((np.array(rates) * (1. - np.array(rates)))/n_trials)
plt.figure()
plt.errorbar(deck_sizes, rates, yerr=errorbars, label="MC analysis, %d trials" % n_trials)
plt.title("Analysis of pair occurance in %d-card hands with %d of each card in deck" % (opening_hand_size, n_card_copies))
plt.ylabel("Fraction of hands containing both cards")
plt.xlabel("Deck size (cards)")
plt.xticks(deck_sizes)

plt.scatter(deck_sizes, gt_rates, c='red', marker=(5, 2), s=80, label="Theoretical expected")
plt.legend()
plt.show()