<jupyter_start><jupyter_text>Setup<jupyter_code>%load_ext autoreload
%autoreload 2

from functools import partial
import math

from sklearn.gaussian_process import GaussianProcessRegressor, GaussianProcessClassifier
from sklearn.model_selection import cross_val_predict
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.neighbors import KernelDensity
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.base import clone
from sklearn.neighbors import KNeighborsRegressor
from xgboost import XGBRegressor

import numpy as np
import seaborn as sns
import scipy
import pandas as pd

from machine_learning.ml_data import MLData
from machine_learning.sklearn import year_cv_split, ColumnDropper, EloRegressor, TeammatchToMatchConverter
from machine_learning.ml_estimators import StackingEstimator
from machine_learning.settings import TEAM_NAMES, ROUND_TYPES, VENUES, CATEGORY_COLS
from machine_learning.ml_estimators.stacking_estimator import ELO_MODEL_COLS

YEAR_RANGE = (2005, 2015)
data = MLData()
data.data
leg_data = MLData(
    pipeline="legacy",
    data_set="legacy_model_data",
)
leg_data.data
se = StackingEstimator()
se
X_train, y_train = data.train_data
X_test, y_test = data.test_data
def _make_cv_prediction(estimator, X, y, train, test):
    X_train, y_train = X[train], y[train]
    X_test = X[test]
    estimator.fit(X_train, y_train)
    return estimator.predict(X_test)

def time_series_cross_val_predict(estimator, X, y, cv):
    make_cv_prediction = partial(_make_cv_prediction, estimator, X, y)

    return np.array([make_cv_prediction(train, test) for train, test in cv])
LOG_BASE = 2
MIN_VAL = 1 * 10**-10

def _calculate_bits(row):
    if row['home_pred'] > row['away_pred']:
        predicted_win_proba = row['home_pred']
        predicted_home_win = True
    else:
        predicted_win_proba = row['away_pred']
        predicted_home_win = False

    correct = (predicted_home_win and row['home_win']) or (not predicted_home_win and not row['home_win'])

    if row['draw']:
        return 1 + (0.5 * math.log(max(predicted_win_proba * (1 - predicted_win_proba), MIN_VAL), LOG_BASE))

    if correct:
        return 1 + math.log(max(row['home_pred'], MIN_VAL), LOG_BASE)

    return 1 + math.log(max(1 - predicted_win_proba, MIN_VAL), LOG_BASE)


def bits_scorer(estimator, X, y, y_pred=None):
    if y_pred is None:
        y_pred = estimator.predict(X)

    team_match_data_frame = X.assign(y_true=y.to_numpy(), y_pred=y_pred)
    home_match_data_frame = team_match_data_frame.query("at_home == 1").sort_index()
    away_match_data_frame = (
        team_match_data_frame.query("at_home == 0")
        .set_index(["oppo_team", "year", "round_number"])
        .rename_axis([None, None, None])
        .sort_index()
    )
    
    bits_data_frame = pd.DataFrame(
        {
            'home_win': home_match_data_frame["y_true"] > away_match_data_frame["y_true"],
            'draw': home_match_data_frame["y_true"] == away_match_data_frame["y_true"],
            'home_pred': home_match_data_frame["y_pred"],
            'away_pred': away_match_data_frame["y_pred"],
        }
    )
    
    return bits_data_frame.apply(_calculate_bits, axis=1)

def neg_bits_eval_func(X, y_true, y_pred):
    team_match_data_frame = X.assign(y_true=y_true.to_numpy(), y_pred=y_pred)
    home_match_data_frame = team_match_data_frame.query("at_home == 1").sort_index()
    away_match_data_frame = (
        team_match_data_frame.query("at_home == 0")
        .set_index(["oppo_team", "year", "round_number"])
        .rename_axis([None, None, None])
        .sort_index()
    )
    
    bits_data_frame = pd.DataFrame(
        {
            'home_win': home_match_data_frame["y_true"] > away_match_data_frame["y_true"],
            'draw': home_match_data_frame["y_true"] == away_match_data_frame["y_true"],
            'home_pred': home_match_data_frame["y_pred"],
            'away_pred': away_match_data_frame["y_pred"],
        }
    )
    
    # Return negative bits, because XGB minimizes evaluation metric by default
    return 'neg_bits', bits_data_frame.apply(_calculate_bits, axis=1).sum() * -1<jupyter_output><empty_output><jupyter_text>Kernel Density Estimation<jupyter_code>cv_y_preds = time_series_cross_val_predict(se, X_train, y_train, year_cv_split(X_train, YEAR_RANGE))

cv_y_preds
sns.scatterplot(
    x=y_train.loc[(slice(None), slice(YEAR_RANGE[0], YEAR_RANGE[1] - 1), slice(None))].to_numpy(),
    y=np.concatenate(cv_y_preds)
)
sns.distplot(np.concatenate(cv_y_preds))
kd = KernelDensity()

kd.fit(np.concatenate(cv_y_preds[:-1]).reshape(-1, 1))
sample = kd.sample(n_samples=1000)

sns.distplot(sample.reshape(-1, 1))
np.concatenate(cv_y_preds).mean(), sample.mean()
np.concatenate(cv_y_preds).std(), sample.std()
confidences = []

for pred in cv_y_preds[-1]:
    confidences.append(scipy.stats.percentileofscore(sample, pred) / 100)
    
conf_ar = np.array(confidences)
conf_ar
X_val = X_train.query('year == 2014')
y_val = y_train.loc[(slice(None), 2014, slice(None))]

bits = bits_scorer(se, X_val, y_val, y_pred=conf_ar)
bits
bits.sum()<jupyter_output><empty_output><jupyter_text>Bootstrapping<jupyter_code>BOOTSTRAP_YEAR_RANGE = (1966, 2015)

def bootstrap_predictions(estimator, X, y, year_range, n_iter=10, one_year=True):
    """Split data by year for cross-validation for time-series data"""

    test_years = np.random.choice(range(*year_range), n_iter)
    preds = []
    
    for year in test_years:
        X_train = X.query('year < @year')
        X_test = X.query('year == @year') if one_year else X.query('year >= @year')
        y_train = y.loc[X_train.index]
        
        boot_estimator = clone(estimator)
        boot_estimator.fit(X_train, y_train)
        preds.append(boot_estimator.predict(X_test))
        
    return np.concatenate(preds)
boot_preds = bootstrap_predictions(se, X_train, y_train, BOOTSTRAP_YEAR_RANGE)
boot_preds
sns.distplot(boot_preds.reshape(-1, 1))
boot_preds.mean(), boot_preds.std()
confidences = []

for pred in cv_y_preds[-1]:
    confidences.append(scipy.stats.percentileofscore(boot_preds, pred) / 100)
    
boot_conf_ar = np.array(confidences)
boot_conf_ar
bits = bits_scorer(se, X_val, y_val, y_pred=boot_conf_ar)
bits
bits.sum()<jupyter_output><empty_output><jupyter_text>Gaussian Process<jupyter_code>gpc_slice = (slice(None), slice(1995, 2013), slice(None))

gpc = make_pipeline(
    ColumnDropper(cols_to_drop=ELO_MODEL_COLS),
    ColumnTransformer(
        [
            (
                "onehotencoder",
                OneHotEncoder(
                    categories=[TEAM_NAMES, TEAM_NAMES, ROUND_TYPES, VENUES],
                    sparse=False,
                    handle_unknown="ignore",
                ),
                CATEGORY_COLS,
            )
        ],
        remainder=StandardScaler(),
    ),
    GaussianProcessClassifier(copy_X_train=False)
)

gpc.fit(X_train.loc[gpc_slice, :], (y_train.loc[gpc_slice] >= 0).astype(int))
gpc_y_pred = gpc.predict_proba(X_test)[:, 1]

bits = bits_scorer(gpc, X_test, (y_test >= 0).astype(int), y_pred=gpc_y_pred)
bits
bits.sum()
gpc_y_pred<jupyter_output><empty_output><jupyter_text>Estimate mean/error<jupyter_code>dist_slice = (slice(None), slice(YEAR_RANGE[0], YEAR_RANGE[1] - 2), slice(None))

error_model = make_pipeline(
    ColumnDropper(cols_to_drop=ELO_MODEL_COLS),
    ColumnTransformer(
        [
            (
                "onehotencoder",
                OneHotEncoder(
                    categories=[TEAM_NAMES, TEAM_NAMES, ROUND_TYPES, VENUES],
                    sparse=False,
                    handle_unknown="ignore",
                ),
                CATEGORY_COLS,
            )
        ],
        remainder=StandardScaler(),
    ),
    KNeighborsRegressor()
)

validation_error = (
    y_train
    .loc[dist_slice]
    .to_numpy() - np.concatenate(cv_y_preds[:-1])
)**2
error_model.fit(X_train.loc[dist_slice, :], validation_error)

pred_mean = cv_y_preds[-1]
pred_std = error_model.predict(X_val)**0.5

pred_mean.mean(), pred_std.mean()
confidences = []

for pred, std in zip(pred_mean, pred_std):
    rv = scipy.stats.norm(loc=pred, scale=std)
    
    if pred > 0:
        confidences.append(1 - rv.cdf(0))
    else:
        confidences.append(rv.cdf(0))
        
err_conf_ar = np.array(confidences)
err_conf_ar
bits = bits_scorer(se, X_val, y_val, y_pred=err_conf_ar)
bits
bits.sum()<jupyter_output><empty_output><jupyter_text>XGBRegressor<jupyter_code>conf_model = make_pipeline(
    ColumnDropper(cols_to_drop=ELO_MODEL_COLS),
    ColumnTransformer(
        [
            (
                "onehotencoder",
                OneHotEncoder(
                    categories=[TEAM_NAMES, TEAM_NAMES, ROUND_TYPES, VENUES],
                    sparse=False,
                    handle_unknown="ignore",
                ),
                CATEGORY_COLS,
            )
        ],
        remainder=StandardScaler(),
    ),
    XGBRegressor(objective='reg:squarederror')
)

conf_model.fit(
    X_train,
    (y_train >= 0).astype(int),
    xgbregressor__eval_metric=neg_bits_eval_func,
)

xgb_y_val = (y_val >= 0).astype(int)
xgb_y_val
bits = bits_scorer(conf_model, X_val, xgb_y_val)
bits
bits.sum()<jupyter_output><empty_output><jupyter_text>ELO<jupyter_code>elo = make_pipeline(TeammatchToMatchConverter(), EloRegressor())
elo.fit(X_train, y_train)
bits = bits_scorer(elo, X_test, y_test)
bits
bits.sum()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1163910404019259, 'fraction_numerical': 0.009838810969227549, 'mean_word_length': 3.7751124437781107, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 38, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '391965', 'license': 'MIT', 'n_tokens_mistral': 3587, 'n_tokens_neox': 3494, 'n_words': 551, 'path': 'notebooks/2020_season/1.0-prediction-confidence-experiments.ipynb', 'repo_name': 'tipresias/augury', 'size': '9554'}