<jupyter_start><jupyter_text>1. Day of Week Analysis<jupyter_code>data = df.groupby('day_of_week').births.mean().reset_index()
font_dict = {'fontsize' : 20}
fig, ax = plt.subplots(figsize = (10, 5))
sns.barplot(
  data = df,
  x = 'day_of_week',
  y = 'births',
  ci = 99
)

ax.set_xlabel('Day of Week', fontdict = font_dict)
ax.set_xticklabels([
  'Monday',
  'Tuesday',
  'Wednesday',
  'Thursday',
  'Friday',
  'Saturday',
  'Sunday'
])
ax.set_ylabel('Mean Births per Day', fontdict = font_dict)
ax.set_title('Births by Day of Week, US 2000–2014', fontdict = font_dict)
plt.xticks(fontsize = 15)
plt.yticks(fontsize = 15)
fig.savefig('outputs/dow.png')

plt.show()
data = df.groupby('month').births.mean().reset_index()

font_dict = {'fontsize' : 20}
fig, ax = plt.subplots(figsize = (10, 5))
sns.barplot(
  data = df,
  x = 'month',
  y = 'births',
  ci = 99
)

ax.set_xlabel('Month', fontdict = font_dict)
ax.set_xticklabels(months)
ax.set_ylabel('Mean Births per Day', fontdict = font_dict)
ax.set_title('Births by Month, US 2000–2014', fontdict = font_dict)
plt.xticks(fontsize = 15, rotation = 90)
plt.yticks(fontsize = 15)
fig.savefig('outputs/month.png')

plt.show()
data = df.groupby('year').births.mean().reset_index()

font_dict = {'fontsize' : 20}
fig, ax = plt.subplots(figsize = (10, 5))
sns.barplot(
  data = df,
  x = 'year',
  y = 'births',
  ci = 99
)

ax.set_xlabel('Year', fontdict = font_dict)
# ax.set_xticklabels(months)
ax.set_ylabel('Mean Births per Day', fontdict = font_dict)
ax.set_title('Births by Year, US 2000–2014', fontdict = font_dict)
plt.xticks(fontsize = 15, rotation = 90)
plt.yticks(fontsize = 15)
fig.savefig('outputs/year.png')

plt.show()<jupyter_output><empty_output><jupyter_text>2. KL–Div Analysis<jupyter_code>from scipy.special import rel_entr as kl_div
# The Overall Distribution
empirical_p = (
  df.groupby(['month', 'date_of_month'])
    .births
    .sum()
    .rename('p_empirical')
)

empirical_p /= empirical_p.sum()

uniform_p = df.groupby(['month', 'date_of_month']).births.count().rename('p_uniform')
uniform_p = uniform_p / uniform_p.sum()

overall_p = pd.DataFrame(empirical_p).join(uniform_p)
# Compute the Probability for Each Year
yearly_p = (
  df.set_index(['month', 'date_of_month', 'year'])
    .groupby('year')
    .births
    .transform(lambda series : series/series.sum())
    .rename('p_empirical')
)

# Compute the Probability if that year were uniform (either 1/355 or 1/356)
uniform_p = 1 / (~yearly_p.unstack('year').isna()).sum().rename('p_uniform')
yearly_p = pd.DataFrame(yearly_p).join(uniform_p)


# Compute the KL-Divergence
overall_kl_div = kl_div(overall_p.p_empirical, overall_p.p_uniform).sum()
yearly_kl_div  = kl_div(yearly_p.p_empirical, yearly_p.p_uniform).groupby('year').sum()
print(f'Yearly  KL–Div : { yearly_kl_div.mean():.3f} ±{ yearly_kl_div.std():.3f}')
print(f'Overall KL-Div : {overall_kl_div.mean():.3f}')
divs = yearly_kl_div
divs['Overall'] = overall_kl_div
divs = divs.rename('kl_div').reset_index()
fig, ax = plt.subplots(figsize = (10, 3))
sns.barplot(
  data = divs.reset_index(),
  y    = 'kl_div',
  x    = 'year',
  ax = ax,
  linewidth=0,
  orient = 'v',
  color = sns.color_palette()[1]
)

ax.set_ylabel('Kullback–Leibler Divergence')
ax.set_xlabel('Year')
ax.set_title('Deviation of the Birth Distribution from Uniform by Year and Overall')
ax.xaxis.grid(False)
ax.yaxis.grid(False)

fig.savefig('outputs/kl_div.png', bbox_inches = 'tight')

plt.show()<jupyter_output><empty_output><jupyter_text>3. Regression Analysis  3.1 Code Holidays<jupyter_code>from sklearn.preprocessing import OneHotEncoder
def code_holiday(year, month, day, dow):
  '''Code holidays
  
  New Years (plus or minus a day)
  Valentines
  Leap Day
  Memorial Day (and the preceding weekend)
  July 4th (and the closest weekend)
  Labor Day (and the preceding weekend)
  Thanksgiving and the surrounding days
  Christmas and the surrounding days
  
  '''
  
  # New Years
  if (month, day) == (1, 1):
    return 'New Years'
  elif (month, day) == (12, 31):
    return 'New Years Eve'
  elif (month, day) == (1, 2):
    return 'Day After New Years'
  
  # Valentines
  if (month, day) == (2, 14):
    return 'Valentines'
  
  # Leap Day
  if (month, day) == (2, 29):
    return 'Leap Day'
  
  # Memorial Day (Last Monday)
  if month == 5:
    if dow == 6   and day + 9 > 31 and day < 30: # Saturday before the last Monday
      return 'Memorial Day Weekend'
    elif dow == 7 and day + 8 > 31 and day < 31: # Sunday before the last Monday
      return 'Memorial Day Weekend'
    elif dow == 1 and day + 7 > 31: # last Monday
      return 'Memorial Day'

  # July 4 and any adjacent weekends
  #  No weekends included if July 4 is Tue–Thur
  if month == 7:
    if day == 4:
      return 'July 4'
    if ((day == 2 and dow == 6) or      
        (day == 3 and dow in [6,7]) or  
        (day == 5 and dow in [6,7]) or
        (day == 6 and dow == 7)):
      return 'July 4 Weekend'
    
  # Labor Day (first Monday of Sep)
  if month == 9:
    if dow == 1 and day - 7 < 1:
      return 'Labor Day'
  if month == 8 and dow in [6, 7] and day + (8-dow) in [32, 33]:
    return 'Labor Day Weekend'
  elif month == 9 and dow in [6, 7] and day - (dow-1) <= 0:
    return 'Labor Day Weekend'
  
  # Thanksgiving (4th Thursday in Nov)
  #  also includes contiguous Wed and Friday and Sat/Sun 
  if month in [11,12]:
    nearest_thur = day + (4 - dow)
    # Mondays projected forwards, Sundays backwards
    if month == 12 and nearest_thur <= 0:
      month = 11
      nearest_thur += 30
      
    if month == 11 and nearest_thur - 21 > 0 and nearest_thur - 21 <= 7:
      if dow == 4:
        return 'Thanksgiving'
      elif dow in [5, 6, 7]:
        return 'Thanksgiving Weekend'
      elif dow == 3:
        return 'Thanksgiving Wednesday'
      
  # Christmas (Dec 25)
  #  also includes: the eve, the day after
  if month == 12:
    if day == 25:
      return 'Christmas'
    elif day == 24:
      return 'Christmas Eve'
    elif day == 26:
      return 'Day After Christmas'
    
  return 'Regular Day'
holidays = df.apply(
  lambda row : code_holiday(row.year, row.month, row.date_of_month, row.day_of_week)
, axis = 'columns')

holiday_encode = OneHotEncoder(drop = ['Regular Day'], sparse = False)
holiday_df = holiday_encode.fit_transform(pd.DataFrame(holidays))
holiday_df = (
  pd.DataFrame(
      holiday_df,
      columns=np.delete(holiday_encode.categories_[0].copy(), holiday_encode.drop_idx_[0]),
      index=df.index
   ).astype(int)
)
print(f'{len(holidays.unique())-1} Holidays')
# # Look at the holidays coded
# pd.options.display.max_rows = 400

# pd.concat([
#   df,
#   holidays
# ], axis = 'columns')[holiday_df.sum(axis = 'columns') == 1]<jupyter_output><empty_output><jupyter_text>3.2 Code the Month and Day of Week<jupyter_code>def encode_col(data, english_names):
  names = dict(enumerate(english_names))
  def name(n):
    return names[n-1]
  
  df = data.apply(name)
  ohe  = OneHotEncoder(drop = 'first', categories = [english_names], sparse = False)
  df = ohe.fit_transform(pd.DataFrame(df))
    
  return pd.DataFrame(
      df,
      columns=np.delete(ohe.categories_[0].copy(), ohe.drop_idx_[0]),
      index=data.index
   ).astype(int)

month_df = encode_col(df.month, months)
dow_df   = encode_col(df.day_of_week, days)<jupyter_output><empty_output><jupyter_text>3.3 Fit the Model<jupyter_code>import statsmodels.api as sm
exog = pd.concat([
  df.year - 2000,
  month_df,
  dow_df,
  holiday_df
], axis = 'columns')

exog = sm.add_constant(exog)
exog.rename({
  'const' : 'Intercept',
  'year'  : 'Year'
},  axis = 'columns', inplace = True)

endog = df.births
model = sm.OLS(endog, exog)
result = model.fit()
print(result.summary())
def barplot(coef_names, ax, name, base_level = None):

  coefs = result.params
  coefs = coefs[coefs.index.isin(coef_names)]
  
  ci    = result.conf_int()
  ci    = ci[ci.index.isin(coef_names)]
  
  if base_level != None:
    coefs [base_level] = 0
    ci.loc[base_level] = 0
    
  sns.barplot(
    data = coefs.rename('Coefficient').reset_index(),
    y = 'index',
    x = 'Coefficient',
    ax = ax,
    linewidth = 0,
    color = sns.color_palette()[0]
  )
  
  for y, interval in enumerate(ci.itertuples()):
    ax.plot((interval._1,interval._2),(y,y),'-',color='orange')
    
  ax.set_xlabel('Change in Births per Day')
  ax.set_ylabel('')
  ax.set_title(name)

fig, axes = plt.subplots(
  ncols = 3, 
  figsize = (20,5)
)
for covariates, name, ax in zip(
    [holiday_df.columns, month_df.columns, dow_df.columns],
    ['Holidays', 'Months', 'Day of Week'],
    axes
  ):
  barplot(covariates, ax, name)

fig.savefig('outputs/model_coefs.png', bbox_inches='tight')
plt.show()
resids = pd.concat([
  df,
  result.resid.rename('residual')
], axis = 'columns')
resids = resids.groupby(['month', 'date_of_month']).residual.mean().unstack(-1)
resids = round(resids, -1)
fig, ax = plt.subplots(figsize = (31, 12))
sns.heatmap(resids,
            ax = ax, 
            center = 0, 
            annot = True,
            fmt = '.0f',
            square = True,
            cmap = 'coolwarm',
            cbar = False,
            annot_kws = {'fontdict' : {'fontsize' : 15}}
           )

font_dict = {'fontsize' : 20}
ax.set_facecolor('white')
ax.set_xlabel('Day of Month', font_dict)
ax.set_ylabel('Month', font_dict)
ax.set_yticklabels([
  'January',
  'February',
  'March',
  'April',
  'May',
  'June',
  'July',
  'August',
  'September',
  'October',
  'November',
  'December'
], font_dict, rotation = 0)
ax.set_title('Unexplained Variance in Birthdays', fontdict = {'fontsize' : 40})
plt.xticks(fontsize = 20)
plt.yticks(rotation = 0)
fig.savefig('outputs/resids.png', transparent = True)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14945346817856778, 'fraction_numerical': 0.023495760547553375, 'mean_word_length': 3.5450324976787373, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '338180', 'license': 'MIT', 'n_tokens_mistral': 3960, 'n_tokens_neox': 3679, 'n_words': 897, 'path': 'bonus_material.ipynb', 'repo_name': 'rcharan/birthdays', 'size': '9789'}