monte_carlo_for_NE_of_running.py

"""
Quick script to run a Monte Carlo simulation of the effect of utility assumptions for a blog post about the last play in the super bowl
"""
from __future__ import division

class NashEquilibrium:
    """
    Class for the Nash Equilibrium (just something to hold data)
    """
    def __init__(self, A, B, C, D):
        self.A = A
        self.B = B
        self.C = C
        self.D = D
        self.valid = (A < B < 100) and (100 > C > D) and (0 < A < C) and (B > D > 0)
        self.y = (D- B) / (D - B + A - C)
        self.alpha = ((self.y) * (A-B)+B) / (A)

if __name__ == '__main__':
    import matplotlib.pyplot as plt  # For plots
    import random  # For random sampling
    import scipy  # For summary statistics
    iterations = 10 ** 6
    As = []
    Bs = []
    Cs = []
    Ds = []
    alphas = []
    ys = []
    while len(alphas) < iterations:
        A = random.normalvariate(60, 10)  # Run vs run D
        B = random.normalvariate(95, 5)  # Run vs pass D
        C = random.normalvariate(85, 15)  # Pass vs run D
        D = random.normalvariate(50, 20)  # Pass vs pass D
        ne = NashEquilibrium(A, B, C, D)
        if ne.valid:
            As.append(A)
            Bs.append(B)
            Cs.append(C)
            Ds.append(D)
            alphas.append(ne.alpha)
            ys.append(ne.y)

    plt.figure()
    parameters = {'$A$':As, '$B$': Bs, '$C$': Cs, '$D$': Ds}
    colors = {'$A$': 'blue', '$B$': 'red', '$C$': 'green', '$D$': 'orange'}
    for p in sorted(parameters.keys()):
        plt.hist(parameters[p], alpha=0.25, label=p, color=colors[p], bins=20)
    plt.xlabel('% Success')
    plt.ylabel('Frequency')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
                       ncol=4, mode="expand", borderaxespad=0.)
    plt.savefig('parameters.png', bbox_inches='tight')

    plt.figure()
    plt.hist(alphas, bins=40)
    plt.xlabel('$\\alpha$')
    plt.ylabel('Frequency')
    plt.savefig('alpha.png')

    print 50 * "="
    print 'Mean alpha: {}'.format(scipy.mean(alphas))
    print 'Std alpha: {}'.format(scipy.std(alphas))
    print 'Max alpha: {}'.format(max(alphas))
    print 'Min alpha: {}'.format(min(alphas))
    print 50 * "="

    plt.figure()
    plt.hist(ys, bins=20)
    plt.xlabel('$y$')
    plt.ylabel('Frequency')
    plt.savefig('y.png')

    print 50 * "="
    print 'Mean y: {}'.format(scipy.mean(ys))
    print 'Std y: {}'.format(scipy.std(ys))
    print 'Max y: {}'.format(max(ys))
    print 'Min y: {}'.format(min(ys))
    print 50 * "="