<jupyter_start><jupyter_text>Introduction  MotivationThis notebook follows up `model_options.ipynb`.The key difference is that we filter using the category distance metric (see `bin/wp-get-links` for details), rather than relying solely on the regression to pick relevant articles. Thus, we want to decide what an appropriate category distance threshold is.Our hope is that by adding this filter, we can now finalize the regression algorithm selection and configuration.  Summary(You may want to read this section last, as it refers to the full analysis below.) Q1. Which algorithm?In `model_options`, we boiled down the selection to three choices:1. Lasso, normalized, positive, auto 𝛼 (LA_NPA)2. Elastic net, normalized, positive, auto 𝛼, auto 𝜌 (EN_NPAA)3. Elastic net, normalized, positive, auto 𝛼, manual 𝜌 = ½ (EN_NPAM)The results below suggest three conclusions:1. LA_NPA vs. EN_NPAA:   * EN_NPAA has (probably insignificantly) better RMSE.   * Forecasts look almost identical.   * EN_NPAA chooses a more reasonable-feeling number of articles.   * EN_NPAA is more principled (lasso vs. elastic net).2. EN_NPAM vs. EN_NPAA:   * EN_NPAA has better RMSE.   * Forecasts look almost identical, except EN_NPAM has some spikes in the 2014–2015 season, which probably accounts for the RMSE difference.   * EN_NPAA chooses fewer articles, though EN_NPAM does not feel excessive.   * EN_NPAA is more principles (manual 𝜌 vs. auto).   On balance, **EN_NPAA seems the best choice**, based on principles and article quantity rather than results, which are nearly the same across the board. Q2. What distance threshold?Observations for EN_NPAA at distance threshold 1, 2, 3:* d = 2 is where RMSE reaches its minimum, and it stays more or less the same all the way through d = 8.* d = 2 and 3 have nearly identical-looking predictions.* d = 2 and 3 have very similar articles and coefficient ranking. Of the 10 and 13 articles respectively, 9 are shared and in almost the same order.These suggests that the actual models for d = 2..8 are very similar. Further:* d = 2 or 3 does not have the spikes in 3rd season that d = 1 does. This suggests that the larger number of articles gives a more robust model.* d = 2 or 3 matches the overall shape of the outbreak better than d = 1, though the latter gets the peak intensity more correct in the 1st season.Finally, examining the article counts in `COUNTS` and `COUNTS_CUM`, d = 2 would give very small input sets in some of the sparser cases, while d = 3 seems safer (e.g., "es+Infecciones por clamidias" and "he+שעלת"). On the other hand, Arabic seems to have a shallower category structure, and d = 3 would capture most articles.On balance, **d = 3 seems the better choice**. It performs as well as d = 2, without catching irrelevant articles, and d = 2 seems too few articles in several cases. The Arabic situation is a bit of an unknown, as none of us speak Arabic, but erring on the side of too many articles seems less risky than clearly too few. Q3. What value of 𝜌?In both this notebook and `model_options`, every auto-selected 𝜌 has been 0.9, i.e., mostly lasso. Thus, we will **fix 𝜌 = 0.9** for performance reasons. ConclusionWe select **EN_NPAM with 𝜌 = 0.9**.  Preamble  Imports<jupyter_code>%matplotlib inline

import collections
import gzip
import pickle
import os
import urllib.parse

import numpy as np
import matplotlib as plt
import pandas as pd
import sklearn as sk
import sklearn.linear_model

DATA_PATH = os.environ['WEIRD_AL_YANKOVIC']
plt.rcParams['figure.figsize'] = (12, 4)<jupyter_output><empty_output><jupyter_text>Load, preprocess, and clean data Load and preprocess the truth spreadsheet.<jupyter_code>truth = pd.read_excel(DATA_PATH + '/truth.xlsx', index_col=0)
TRUTH_FLU = truth.loc[:,'us+influenza']  # pull Series
TRUTH_FLU.index = TRUTH_FLU.index.to_period('W-SAT')
TRUTH_FLU.head()<jupyter_output><empty_output><jupyter_text>Load the Wikipedia link data. We convert percent-encoded URLs to Unicode strings for convenience of display.<jupyter_code>def unquote(url):
   (lang, url) = url.split('+', 1)
   url = urllib.parse.unquote(url)
   url = url.replace('_', ' ')
   return (lang + '+' + url)
raw_graph = pickle.load(gzip.open(DATA_PATH + '/articles/wiki-graph.pkl.gz'))
GRAPH = dict()
for root in raw_graph.keys():
   unroot = unquote(root)
   GRAPH[unroot] = { unquote(a): d for (a, d) in raw_graph[root].items() }<jupyter_output><empty_output><jupyter_text>Load all the time series. Most of the 4,299 identified articles were in the data set.Note that in contrast to `model_options`, we do not remove any time series by the fraction that they are zero. The results seem good anyway. This filter also may not apply well to the main experiment, because the training periods are often not long enough to make it meaningful.<jupyter_code>TS_ALL = pd.read_csv(DATA_PATH + '/tsv/forecasting_W-SAT.norm.tsv',
                     sep='\t', index_col=0, parse_dates=True)
TS_ALL.index = TS_ALL.index.to_period('W-SAT')
TS_ALL.rename(columns=lambda x: unquote(x[:-5]), inplace=True)
len(TS_ALL.columns)
(TS_ALL, TRUTH_FLU) = TS_ALL.align(TRUTH_FLU, axis=0, join='inner')
TRUTH_FLU.plot()
TS_ALL.iloc[:,:5].plot()<jupyter_output><empty_output><jupyter_text>Summarize distance from root Number of articles by distance from each root.<jupyter_code>COUNTS = pd.DataFrame(columns=range(1,9), index=sorted(GRAPH.keys()))
COUNTS.fillna(0, inplace=True)
for (root, leaves) in GRAPH.items():
   for (leaf, dist) in leaves.items():
      COUNTS[dist][root] += 1
COUNTS<jupyter_output><empty_output><jupyter_text>Number of articles of at most the given distance.<jupyter_code>COUNTS_CUM = COUNTS.cumsum(axis=1)
COUNTS_CUM<jupyter_output><empty_output><jupyter_text>Parameter sweep Return the set of articles with maximum category distance from a given root.<jupyter_code>def articles_dist(root, dist):
   return { a for (a, d) in GRAPH[root].items() if d <= dist }<jupyter_output><empty_output><jupyter_text>Return time series for articles with a maximum category distance from the given root.<jupyter_code>def select_by_distance(root, d):
   keep_cols = articles_dist(root, d)
   return TS_ALL.filter(items=keep_cols, axis=1)
select_by_distance('en+Influenza', 1).head()<jupyter_output><empty_output><jupyter_text>Fit function. The core is the same as the `model_options` one, with non-constant training series set and a richer summary.<jupyter_code>def fit(root, train_week_ct, d, alg, plot=True):
   ts_all = select_by_distance(root, d)
   ts_train = ts_all.iloc[:train_week_ct,:]
   truth_train = TRUTH_FLU.iloc[:train_week_ct]
   m = alg.fit(ts_train, truth_train)
   m.input_ct = len(ts_all.columns)
   pred = m.predict(ts_all)
   pred_s = pd.Series(pred, index=TRUTH_FLU.index)
   m.r = TRUTH_FLU.corr(pred_s)
   m.rmse = ((TRUTH_FLU - pred_s)**2).mean()
   m.nonzero = np.count_nonzero(m.coef_)
   if (not hasattr(m, 'l1_ratio_')):
      m.l1_ratio_ = -1
   # this is just a line to show how long the training period is
   train_period = TRUTH_FLU.iloc[:train_week_ct].copy(True)
   train_period[:] = 0
   if (plot):
      pd.DataFrame({'truth':       TRUTH_FLU,
                    'prediction':  pred,
                    'training pd': train_period}).plot(ylim=(-1,9))
   sumry = pd.DataFrame({'coefs': m.coef_,
                         'coefs_abs': np.abs(m.coef_)},
                        index=ts_all.columns)
   sumry.sort_values(by='coefs_abs', ascending=False, inplace=True)
   sumry = sumry.loc[:, 'coefs']
   for a in ('intercept_', 'alpha_', 'l1_ratio_', 'nonzero', 'rmse', 'r', 'input_ct'):
      try:
         sumry = pd.Series([getattr(m, a)], index=[a]).append(sumry)
      except AttributeError:
         pass
   return (m, pred, sumry)<jupyter_output><empty_output><jupyter_text>Which 𝛼 and 𝜌 to explore? Same as `model_options`.<jupyter_code>ALPHAS = np.logspace(-15, 2, 25)
RHOS = np.linspace(0.1, 0.9, 9)<jupyter_output><empty_output><jupyter_text>Try all distance filters and summarize the result in a table.<jupyter_code>def fit_summary(root, label, train_week_ct, alg, **kwargs):
   result = pd.DataFrame(columns=[[label] * 4,
                                  ['input_ct', 'rmse', 'rho', 'nonzero']],
                         index=range(1, 9))
   preds = dict()
   for d in range(1, 9):
      (m, preds[d], sumry) = fit(root, train_week_ct, d, alg(**kwargs), plot=False)
      result.loc[d,:] = (m.input_ct, m.rmse, m.l1_ratio_, m.nonzero)
   return (result, preds)<jupyter_output><empty_output><jupyter_text>Lasso, normalized, positive, auto 𝛼<jupyter_code>la_npa = fit_summary('en+Influenza', 'la_npa', 104, sk.linear_model.LassoCV,
                      normalize=True, positive=True, alphas=ALPHAS,
                      max_iter=1e5, selection='random', n_jobs=-1)
la_npa[0]
(m, _, s) = fit('en+Influenza', 104, 1,
                sk.linear_model.LassoCV(normalize=True, positive=True, alphas=ALPHAS,
                                        max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)
(m, _, s) = fit('en+Influenza', 104, 2,
                sk.linear_model.LassoCV(normalize=True, positive=True, alphas=ALPHAS,
                                        max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)
(m, _, s) = fit('en+Influenza', 104, 3,
                sk.linear_model.LassoCV(normalize=True, positive=True, alphas=ALPHAS,
                                        max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)<jupyter_output><empty_output><jupyter_text>Elastic net, normalized, positive, auto 𝛼, auto 𝜌<jupyter_code>en_npaa = fit_summary('en+Influenza', 'en_npaa', 104, sk.linear_model.ElasticNetCV,
                      normalize=True, positive=True, alphas=ALPHAS, l1_ratio=RHOS,
                      max_iter=1e5, selection='random', n_jobs=-1)
en_npaa[0]
(m, _, s) = fit('en+Influenza', 104, 1,
                sk.linear_model.ElasticNetCV(normalize=True, positive=True,
                                             alphas=ALPHAS, l1_ratio=RHOS,
                                             max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)
(m, _, s) = fit('en+Influenza', 104, 2,
                sk.linear_model.ElasticNetCV(normalize=True, positive=True,
                                             alphas=ALPHAS, l1_ratio=RHOS,
                                             max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)
(m, _, s) = fit('en+Influenza', 104, 3,
                sk.linear_model.ElasticNetCV(normalize=True, positive=True,
                                             alphas=ALPHAS, l1_ratio=RHOS,
                                             max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)<jupyter_output><empty_output><jupyter_text>Elastic net, normalized, positive, auto 𝛼, manual 𝜌 = ½<jupyter_code>en_npam = fit_summary('en+Influenza', 'en_npam', 104, sk.linear_model.ElasticNetCV,
                      normalize=True, positive=True, alphas=ALPHAS, l1_ratio=0.5,
                      max_iter=1e5, selection='random', n_jobs=-1)
en_npam[0]
(m, _, s) = fit('en+Influenza', 104, 1,
                sk.linear_model.ElasticNetCV(normalize=True, positive=True,
                                             alphas=ALPHAS, l1_ratio=0.5,
                                             max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)
(m, _, s) = fit('en+Influenza', 104, 2,
                sk.linear_model.ElasticNetCV(normalize=True, positive=True,
                                             alphas=ALPHAS, l1_ratio=0.5,
                                             max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)
(m, _, s) = fit('en+Influenza', 104, 3,
                sk.linear_model.ElasticNetCV(normalize=True, positive=True,
                                             alphas=ALPHAS, l1_ratio=0.5,
                                             max_iter=1e5, selection='random', n_jobs=-1))
s.head(27)<jupyter_output><empty_output><jupyter_text>Summary All the result tables next to one another.<jupyter_code>pd.concat([la_npa[0], en_npaa[0], en_npam[0]], axis=1)<jupyter_output><empty_output><jupyter_text>Plot the predictions by distance filter next to one another.<jupyter_code>def plot(data, ds):
   for d in ds:
      D = collections.OrderedDict([('truth', TRUTH_FLU)])
      D[d] = data[1][d]
      pd.DataFrame(D).plot(figsize=(12,3))
plot(la_npa, range(1, 4))
plot(en_npaa, range(1, 4))
plot(en_npam, range(1, 4))<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup