goalpermin = 18/70.
goalpermin

"""
Explanation: There is one missing team in 2014: Team 17 Flesh Eating Virus.  They went 5-2 with a total plus/minus of +11.  So a per-game +/- of +1.57.  No need to worry about this missing data point.
Part 2: Convert average plus/minus into likelihood of winning.
If I know the average number of goals scored in a game and I know the length of each game, then I can get an estimate of the average scoring rate.
Suppose average number of goals scored per game is 18. And that game length is 70 minutes.
End of explanation
"""


ok = [poisson(goalpermin, 70).sum() for i in range(171)]

sns.distplot(ok, kde=False)

"""
Explanation: Then the average number of goals per minute is 0.257 goals/minute.
End of explanation
"""


18 / 70.

0.257 / 2

"""
Explanation: That looks like a pretty reasonable distribution of point totals per game.
Next, the question we are really interested in: suppose you are on a team that has a plus/minus average per game of -5.  If the halves are split evenly in time, then at halftime you are typically behind by 2.5 points, and there are another 18 or so points to be played.  
I am missing something here.  Need some measure of the variance of plus/minus values.  Time to sleep on it.
End of explanation
"""


gpm1 = 18 / 2 - 3.5

gpm2 = 18/2 + 1.5

gpm1 / 70.

gpm2/70.

ok1 = [poisson(gpm1) for i in range(171)]

ok2 = [poisson(gpm2) for i in range(171)]

sns.distplot(ok1, kde=False, bins=range(25))
sns.distplot(ok2, kde=False, bins=range(25))

okdiff = np.array(ok1) - np.array(ok2)
oksum = np.array(ok1) + np.array(ok2)

sns.distplot(okdiff, kde=False)

sns.distplot(oksum)

np.mean(okdiff)

ok = [poisson(gpm1, 70).sum() + poisson(gpm2, 70).sum() for i in range(171)]

np.mean(ok)

def underdogwin(gpm1, gpm2, remaining_time, thresh):
    ok1 = [poisson(gpm1, 70 - remaining_time).sum() for i in range(100)]
    ok2 = [poisson(gpm2, 70 - remaining_time).sum() for i in range(100)]
    okdiff = np.array(ok1) - np.array(ok2)
    over5 = okdiff > thresh
    return len(okdiff[over5])

nsim = 100
sim_come0 = []
sim_come1 = []
sim_come2 = []
sim_come3 = []
for isim in tqdm(range(nsim)):
    come_from_behind0 = []
    come_from_behind1 = []
    come_from_behind2 = []
    come_from_behind3 = []
    remaining_times = range(70)
    for remaining_time in remaining_times:
        gpm1 = 18 / 2 / 70.
        gpm2 = 18 / 2 / 70.
        thresh = 0
        wins = underdogwin(gpm1, gpm2, remaining_time, thresh)
        come_from_behind0.append(wins)
        
        thresh = 1
        wins = underdogwin(gpm1, gpm2, remaining_time, thresh)
        come_from_behind1.append(wins)
        
        thresh = 2
        wins = underdogwin(gpm1, gpm2, remaining_time, thresh)
        come_from_behind2.append(wins)
        
        thresh = 3
        wins = underdogwin(gpm1, gpm2, remaining_time, thresh)
        come_from_behind3.append(wins)
    sim_come0.append(come_from_behind0)
    sim_come1.append(come_from_behind1)
    sim_come2.append(come_from_behind2)
    sim_come3.append(come_from_behind3)
come_mean0 = np.mean(sim_come0, axis=0)
come_mean1 = np.mean(sim_come1, axis=0)
come_mean2 = np.mean(sim_come2, axis=0)
come_mean3 = np.mean(sim_come3, axis=0)
come_std0 = np.std(sim_come0, axis=0)
come_std1 = np.std(sim_come1, axis=0)
come_std2 = np.std(sim_come2, axis=0)
come_std3 = np.std(sim_come3, axis=0)

nsim = 100
sim_comes = []
threshes = range(4)
for thresh in threshes:
    sim_come = []
    for isim in tqdm(range(nsim)):
        come_from_behind = []
        remaining_times = range(70)
        for remaining_time in remaining_times:
            gpm1 = 18 / 2 / 70.
            gpm2 = 18 / 2 / 70.
            wins = underdogwin(gpm1, gpm2, remaining_time, thresh)
            come_from_behind.append(wins)
        sim_come.append(come_from_behind)
    sim_comes.append(sim_come)
come_mean0 = np.mean(sim_comes[0], axis=0)
come_mean1 = np.mean(sim_comes[1], axis=0)
come_mean2 = np.mean(sim_comes[2], axis=0)
come_mean3 = np.mean(sim_comes[3], axis=0)
come_std0 = np.std(sim_comes[0], axis=0)
come_std1 = np.std(sim_comes[1], axis=0)
come_std2 = np.std(sim_comes[2], axis=0)
come_std3 = np.std(sim_comes[3], axis=0)

nsim = 100
sim_comes = []
threshes = range(4)
for thresh in threshes:
    sim_come = []
    for isim in tqdm(range(nsim)):
        come_from_behind = []
        remaining_times = range(70)
        for remaining_time in remaining_times:
            gpm1 = 18 / 2 / 70. - 3.5 / 70
            gpm2 = 18 / 2 / 70. + 1.5 / 70
            wins = underdogwin(gpm1, gpm2, remaining_time, thresh)
            come_from_behind.append(wins)
        sim_come.append(come_from_behind)
    sim_comes.append(sim_come)
come_mean0_bad = np.mean(sim_comes[0], axis=0)
come_mean1_bad = np.mean(sim_comes[1], axis=0)
come_mean2_bad = np.mean(sim_comes[2], axis=0)
come_mean3_bad = np.mean(sim_comes[3], axis=0)
come_std0_bad = np.std(sim_comes[0], axis=0)
come_std1_bad = np.std(sim_comes[1], axis=0)
come_std2_bad = np.std(sim_comes[2], axis=0)
come_std3_bad = np.std(sim_comes[3], axis=0)

sns.set_context('poster')
sns.set_style('white')
f, ax = plt.subplots(figsize=(6, 5))
xarr = np.arange(len(come_mean2_bad), 0, -1)
yarr = gaussian_filter1d(come_mean2, 3)
sigarr = gaussian_filter1d(come_std2, 3)
y1 = yarr - sigarr
y2 = yarr + sigarr
# plt.fill_between(xarr, y1, y2, color='gray', alpha=0.4)
plt.plot(xarr, yarr, label='Team A equal to Team B')


yarr = gaussian_filter1d(come_mean2_bad, 3)
sigarr = gaussian_filter1d(come_std2_bad, 3)
y1 = yarr - sigarr
y2 = yarr + sigarr
# plt.fill_between(xarr, y1, y2, color='gray', alpha=0.4)
plt.plot(xarr, yarr, label='Team A much worse than Team B')

plt.ylim([0, 100])
plt.xlim([70,0])
plt.xlabel('Time Remaining [minutes]')
plt.ylabel('Percent Chance Team A Wins By Three')
plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(figures_dir, 'WinBy3Probability'))
# plt.plot(come_mean_bad)

sns.set_context('poster')
# xarr = range(len(come_mean2))
# yarr = gaussian_filter1d(come_mean2, 3)
# sigarr = gaussian_filter1d(come_std2, 3)
# y1 = yarr - sigarr
# y2 = yarr + sigarr
# plt.fill_between(xarr, y1, y2, color='gray', alpha=0.5)
# plt.plot(xarr, yarr)

xarr = range(len(come_mean0_bad))
yarr = gaussian_filter1d(come_mean0_bad, 3)
sigarr = gaussian_filter1d(come_std0_bad, 3)
y1 = yarr - sigarr
y2 = yarr + sigarr
plt.fill_between(xarr, y1, y2, color='gray', alpha=0.5)
plt.plot(xarr, yarr)
# plt.plot(come_mean_bad)

sph_index = (ratings.allteams['season'] == 'Spring') & \
            (ratings.allteams['type'] == 'Hat') & \
            (ratings.allteams['divname'] == 'JP Mixed (4/3)') & \
            (ratings.allteams['year'] >= 2010)
sph = ratings.allteams[sph_index]

sns.distplot(sph['n_exp_rating'], kde=False, bins=10)
sns.distplot(sph['n_cap_rating'], kde=False, bins=10)
sns.distplot(sph['n_capexp_rating'], kde=False, bins=10)

sph['n_exp_rating'].mean()

sph['n_cap_rating'].mean()

sph['n_capexp_rating'].mean()

sph['n_exp_rating'].median()

sph['n_cap_rating'].median()

sph['n_capexp_rating'].median()

.75*16

.69*16

.84*16

"""
Explanation: An average team playing against an average team has an expected goal scoring rate of 0.257 / 2 = 0.128 goals/minute. A team that averages -5 plus/minus per game gives up 2.5 goals more pre 70 minutes and scores 2.5 goals less per 70 minutes while playing an average team.
End of explanation
"""


goalpermin

avgoff = []
for isim in range(171):
    ok1 = [poisson(goalpermin/2, 70).sum() for i in range(7)]
    ok2 = [poisson(goalpermin/2, 70).sum() for i in range(7)]
    off = np.array(ok1) - np.array(ok2)
    avgoff.append(off.mean())

sns.set_context('poster')
sns.set_style('white')
fig, ax = plt.subplots(figsize=(6, 5))
sns.distplot(avgoff, kde=False, bins=range(-10, 11), color='gray', label='Equal Skill')
sns.distplot(sph['plusminus'], kde=False, bins=range(-10,11), label='Observed')
plt.ylabel('Number of Teams')
plt.xlabel('Average Plus/Minus per Game')
plt.tight_layout()
plt.legend()
plt.savefig(os.path.join(figures_dir, 'PlusMinusDistribution'))

"""
Explanation: Simulation of two evenly matched teams.  What is the distribution of 171 instances of observed average point differential over 7 games?
End of explanation
"""


wins = []
for isim in range(171):
    ok1 = [poisson(goalpermin/2, 70).sum() for i in range(7)]
    ok2 = [poisson(goalpermin/2, 70).sum() for i in range(7)]
    off = np.array(ok1) - np.array(ok2)
    win_index = off > 0
    wins.append(len(off[win_index]))

sns.distplot(wins, kde=False, bins=range(8))

np.array(wins).sum()

7 * 171

564 / 1197.

1/7.

n1 = len(np.array(wins)[np.array(wins) < 1])

n1

2 / 171.

sns.set_context('poster')
sns.set_style('white')
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
sns.distplot(avgoff, kde=False, bins=range(-10, 11), color='gray', label='Equal Skill', ax=ax)
sns.distplot(sph['plusminus'], kde=False, bins=range(-10,11), label='Observed', ax=ax)
ax.set_ylabel('Number of Teams')
ax.set_xlabel('Average Plus/Minus per Game')
ax.legend()

ax = axes[1]
xarr = np.arange(len(come_mean2_bad), 0, -1)
yarr = gaussian_filter1d(come_mean2, 3)
sigarr = gaussian_filter1d(come_std2, 3)
y1 = yarr - sigarr
y2 = yarr + sigarr
# plt.fill_between(xarr, y1, y2, color='gray', alpha=0.4)
ax.plot(xarr, yarr, label='Team A equal to Team B')


yarr = gaussian_filter1d(come_mean2_bad, 3)
sigarr = gaussian_filter1d(come_std2_bad, 3)
y1 = yarr - sigarr
y2 = yarr + sigarr
# plt.fill_between(xarr, y1, y2, color='gray', alpha=0.4)
ax.plot(xarr, yarr, label='Team A much worse than Team B')

ax.set_ylim([0, 100])
ax.set_xlim([70,0])
ax.set_xlabel('Time Remaining [minutes]')
ax.set_ylabel('Percent Chance Team A Wins By Three')
ax.legend()
plt.tight_layout(w_pad=2)
plt.savefig(os.path.join(figures_dir, 'PlusMinusDistribution_WinBy3Probability'))
# plt.plot(come_mean_bad)

"""
Explanation: Given equal teams, what is the likelihood of winning fewer than 10% of your games?
End of explanation
"""


ways_to_win_one_or_less = 7 + 1

total_possible_outcomes = 2. ** 7

probability = ways_to_win_one_or_less / total_possible_outcomes

probability

1 / 7.

7 * 5 * 4 / 3 / 2. / 2**7

corepart = 1.396 + 1.476 + 1.484 + 1.363 + 1.082 + 1.038

corepart / 6.

lcp = 786. + 800 + 819 + 899 + 1032 + 947
lcp / 6

rcp = 342. + 306 + 361 + 427 + 440 + 431
rcp / 6

(13651 - 8467) / 8467.

8467 * 1.6

7/11.

leaguelist = {'alabama': 3, 'arkansas': 3, 'alaska': 3, 'arizona': 3, 'california': 10, 'colorado': 7, 'connecticut': 4,  
 'delaware': 1, 'florida': 9, 'georgia': 5, 'hawaii': 1, 'idaho': 1, 'illinois': 4, 'indiana': 8, 'iowa': 5, 
 'kansas': 4, 'kentucky': 2, 'louisiana': 1, 'maine': 3, 'maryland': 5, 'massachusetts': 2, 
 'michigan': 6, 'minnesota': 2, 'mississippi': 2, 'missouri': 2, 'montana': 2, 'nebraska': 2, 'nevada': 2, 
 'new hampshire': 1, 'new jersey': 5, 'new mexico': 3, 'new york': 11, 'north carolina': 8, 'north dakota': 0,
 'ohio': 4, 'oklahoma': 1, 'oregon': 6, 'pennsylvania': 6, 'rhode island': 1, 'south carolina': 5, 'south dakota': 1,
 'tennessee': 4, 'texas': 6, 'utah': 4, 'vermont': 2, 'virginia': 6, 
 'washington': 3, 'west virginia': 0, 'wisconsin': 3, 'wyoming': 0}

totleagues = 0
for key in leaguelist:
    totleagues += leaguelist[key]

totleagues

yarr

200/60.

1/3. * (1800 + 1580 + 1400)

16/6.

.84*16

import seaborn as sns

current_palette = sns.color_palette()
sns.palplot(current_palette)

current_palette

"""
Explanation: An alternative way of investigating the probability of winning one game out of 7 or less: simple probability.
End of explanation
"""


