# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Python Practice Lecture 16 MATH 342W Queens College - K-fold Cross Validation
# ## Author: <NAME>
# ## Date: April 5, 2022
#
# ## Agenda:
# * The K Tradeoff
# * Reducing variance with Cross Validation (i.e. K-fold CV)
#
# ## The K Tradeoff
#
# $K$ determines how large the training set is relative to the test set when you're doing honest validation for an algorithm. For now, let's not use K-fold CV, but only examine one split at a time. Consider this simulated dataset with 50 observations:

# +
# importing dependencies
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error # calculates RMSE
import pandas as pd
import random

# ignore warnings
import warnings
warnings.filterwarnings("ignore")

# Data viz
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_palette(sns.color_palette("colorblind")) # setting color palette
sns.set(rc={"figure.figsize":(10, 6)}) #width=10, #height=6

# +
# setting seed
np.random.seed(2022)

# constants
n = 50
xmin = 0
xmax = 4

# data vecs
x = np.random.uniform(low = xmin, high = xmax, size = n)
y = [2 + 3 * i ** 2 + np.random.normal(0, 0.8) for i in x]

# convert to df
Xy = pd.DataFrame({'y': y, 'x': x})

# snapshot
Xy.head()

# +
from plotnine import ggplot, aes, geom_point

# plot the data
data_plot = ggplot(Xy) + aes(x = x, y = y) + geom_point()
data_plot
# -

# Note how $f(x)$ is quadratic and there is random noise which is "ignorance error". The random noise will be part of generalization error and can never go away.
#
# If we use OLS with no derived features, then we can at most get $h*(x)$. Let's see what $h^*(x) = \beta_0 + \beta_1 x$ truly is. To do this, we imagine we see an absolute ton of data and run OLS on it.

# +
# setting seed
np.random.seed(1738)

n_hidden = 1e6

# data vecs
x_hidden = np.random.uniform(low = xmin, high = xmax, size = n)
y_hidden = [2 + 3 * i ** 2 + np.random.normal(0, 0.8) for i in x_hidden]

# convert to dataframes
x_hidden = pd.DataFrame(x_hidden)
y_hidden = pd.DataFrame(y_hidden)


# h_*
h_star_model = LinearRegression().fit(x_hidden, y_hidden)

print(h_star_model.intercept_[0], h_star_model.coef_[0][0])
# -

# The fact that $\beta = [-6~12]^\top$ can actually be solved with calculus: $\int_0^4 ((2 + 3x^2) - (b_0 + b_1 x))^2 dx$ and solve for $b_0$ and $b_1$ explicitly by minimizing.
#
# Plotting that over $\mathbb{D}$ we obtain

# +
from plotnine import geom_abline

data_plot + geom_abline(
    intercept = h_star_model.intercept_[0],
    slope = h_star_model.coef_[0][0],
    color = 'green')
# -

# That is the best we're going to get. However, $g_{final}$ falls far short of it:

# +
# h_*
g_final_model = LinearRegression().fit(Xy[['x']], Xy[['y']])

print(g_final_model.intercept_[0], g_final_model.coef_[0][0])
# -

# The actual standard error of g_final can be estimated by imagining tons of future observations:

# +
y_hat_g_final = g_final_model.predict(x_hidden)

gen_error_true = np.std(y_hidden - y_hat_g_final)
gen_error_true[0]
# -

# The model $g$ can vary quite a bit as we subsample $\mathbb{D}$ which is what happens when you do train-test splits. It varies a lot because there is large misspecification error. If the model was correctly specified, the results of everything that follows will be less impressive. But in the real world - is your model ever correctly specified? is $f \in \mathcal{H}$?? NO. So this is more realistic.
#
# Now let's let K be small. Let K = 2 meaning even 50-50 split of test and train.
#
# Sets, initialized with paranthesis, not square brackets, allow us to reduce one list from the other using the subtract operator.

# +
# setting seed
random.seed(17388) # two different libraries to set seeds for

# constants
K = 2
prop_train = (K - 1) / K
n_train = int(prop_train * n)

# train test split indices
indices = set([i for i in range(n)])
index_train = set(random.sample(indices, n_train))
index_test = indices - index_train

# split based on indices
x_train = pd.DataFrame([x[i] for i in index_train]) # cast from set to list
y_train = pd.DataFrame([y[i] for i in index_train])

x_test = pd.DataFrame([x[i] for i in index_test])
y_test = pd.DataFrame([y[i] for i in index_test])

# fitting model
g_mod = LinearRegression().fit(x_train, y_train)

# yhat via test set
y_hat_g_test = g_mod.predict(x_test)

# yhat via train set
y_hat_g_train = g_mod.predict(x_train)

# standard error
gen_error_K_2_test = np.std(y_test - y_hat_g_test) # test
gen_error_K_2_train = np.std(y_train - y_hat_g_train) # train

print(gen_error_true[0]) # true standard error
print(gen_error_K_2_test[0]) # tested standard error
# -

# Although I cooked the books by setting the seed, this realization makes sense. If K=2, I build the model g with half the data than the model g_final. Less data to train on => higher generalization error. How about if K is large. Let's say $K = \frac{n}{2}$ meaning $n_{train} = 48$ and $n_{test} = 2$.

# +
# setting seed
np.random.seed(7777)
random.seed(7777) # two different libraries to set seeds for

# constants
K = n / 2
prop_train = (K - 1) / K
n_train = int(prop_train * n)

# train test split indices
indices = set([i for i in range(n)])
index_train = set(random.sample(indices, n_train))
index_test = indices - index_train

# split based on indices
x_train = pd.DataFrame([x[i] for i in index_train]) # cast from set to list
y_train = pd.DataFrame([y[i] for i in index_train])

x_test = pd.DataFrame([x[i] for i in index_test])
y_test = pd.DataFrame([y[i] for i in index_test])

# fitting model
g_mod = LinearRegression().fit(x_train, y_train)

# yhat via test set
y_hat_g_test = g_mod.predict(x_test)

# yhat via train set
y_hat_g_train = g_mod.predict(x_train)

# standard error
gen_error_K_2_test = np.std(y_test - y_hat_g_test) # test
gen_error_K_2_train = np.std(y_train - y_hat_g_train) # train

print(gen_error_true[0]) # true standard error
print(gen_error_K_2_test[0]) # tested standard error
print(gen_error_K_2_train[0]) # train standard error
# -

# Although I cooked the books again by setting the seed, this also makes sense. More data to train on means less error but still more error than all the data. In reality, there is massive variance over specific splits! Let's run the simulation with these two K values many times.
#
# While we're at it, let's do all K's! Well, what are all the valid K's? If you want to keep the sizes the same, any factorization of n except the trivial 1 since n = 1 * n. A $K = 1$ would mean there's no split! How to find divisors? Of course a package for this.

# +
# import divisors() method from sympy
from sympy import divisors

Kuniques = list(set(divisors(n)) - set([1]))
Kuniques
# -

# But should we also include the trivial n? Yes K = n is indeed a valid divisor. And this type of CV is called the "leave one out cross validation" (LOOCV). Now we compute the errors over K:

# +
# constants
Nsim_per_K = 2000 # simulations per K fold
num_Kuniques = len(Kuniques)
Ks = Kuniques * Nsim_per_K

results = pd.DataFrame({'standard_error': [0 for i in range(Nsim_per_K * num_Kuniques)],
                        'K': [0 for i in range(Nsim_per_K * num_Kuniques)]})

for i in range(len(Ks)):
    # constants
    K = Ks[i]
    prop_train = (K - 1) / K
    n_train = int(prop_train * n)

    # train test split indices
    indices = set([i for i in range(n)])
    index_train = set(random.sample(indices, n_train))
    index_test = indices - index_train

    # split based on indices
    X_train = pd.DataFrame([x[i] for i in index_train]) # cast from set to list
    y_train = pd.DataFrame([y[i] for i in index_train])

    X_test = pd.DataFrame([x[i] for i in index_test])
    y_test = pd.DataFrame([y[i] for i in index_test])

    # fitting model
    g_mod = LinearRegression().fit(X_train, y_train)
    
    # get predictions
    y_hat_g = g_mod.predict(X_test)
    
    # get error metric
    if len(y_test) == 1:
        g_s_e = np.absolute([float(y_test.iloc[i]) - y_hat_g[i] for i in range(len(y_test))])
    else:
        g_s_e = np.std([float(y_test.iloc[i]) - y_hat_g[i] for i in range(len(y_test))])
    
    # update results
    results.loc[i] = np.asarray([g_s_e, K])
    
# output results
results
# -

# What are the variabilities? Let's take the average error over each simulated split. 

# +
# mean of the standard errors grouped by K
results_summary = pd.DataFrame({'K': list(results.groupby(['K']).mean().index),
                                 'K_avg': list(results.groupby(['K']).mean()['standard_error']),
                                 'K_SE': list(results.groupby(['K']).std()['standard_error'])
                                }).set_index('K')

results_summary

# +
sns.displot(results,
             x = 'standard_error',
             hue = 'K',
             kind = 'kde',
             multiple="stack",
             palette = 'bright',
             alpha = 0.9,
             aspect = 1.5) # width = height * aspect

plt.axvline(label='Mean 2.0', x=results_summary.iloc[0, 0], color='blue')
plt.axvline(label='Mean 5.0', x=results_summary.iloc[1, 0], color='orange')
plt.axvline(label='Mean 10.0', x=results_summary.iloc[2, 0], color='green')
plt.axvline(label='Mean 25.0', x=results_summary.iloc[3, 0], color='red')
plt.axvline(label='Mean 50.0', x=results_summary.iloc[4, 0], color='purple')
# -

# The main takeaways are
#
# 1. The std err of generalization error estimate is much lower for low K than high K
#
# With high K, the test set is small meaning the estimate has high variance; with low K, the test set is large meaning you can measure it with low variance.
#
# 2. The average of generalization error estimate is lower for high K than low K
#
# With high K, the training set is large meaning $g$ is closer to g_final and thus has higher expected accuracy; with low K, the training set is small meaning $g$ is further from g_final and thus has lower expected accuracy.
#
# Thus, the tradeoff is bias vs. variance. There are many similar tradeoffs in statistics. We will see one later when we do machine learning.
#
# Is the estimates' accuracy for what we really care about? No, the generalization error of g_final which we picture below:

# +
sns.displot(results,
             x = 'standard_error',
             hue = 'K',
             kind = 'kde',
             multiple="stack",
             palette = 'bright',
             alpha = 0.9,
             aspect = 1.5) # width = height * aspect

plt.axvline(label='Mean 2.0', x=results_summary.iloc[0, 0], color='blue')
plt.axvline(label='Mean 5.0', x=results_summary.iloc[1, 0], color='orange')
plt.axvline(label='Mean 10.0', x=results_summary.iloc[2, 0], color='green')
plt.axvline(label='Mean 25.0', x=results_summary.iloc[3, 0], color='red')
plt.axvline(label='Mean 50.0', x=results_summary.iloc[4, 0], color='purple')
plt.axvline(label='CV', x=gen_error_true[0], color='white')

# -

# Remember, g_final's error should be lower than both averages since it uses all the data. But we see above it's higher!
#
# So what happened? Simple... we are mixing apples and oranges. We calculated that white line by looking at one million future observations. We calculated the red and blue distributions by looking at our data only which is a random realization of many such datasets! Thus, our generalization errors are biased based on the specific n observations in D we received. We will see that K-fold helps a bit with this. But there is nothing we can do about it beyond that (besides collect more observations). If you get a weird sample, you get a weird sample!
#
# How would we be able to generate the picture we really want to see? We would run this simulation over many datasets and average. That would be a giant simulation. To show that this is the case, go back and change the seed in the first chunk and rerun. You'll see a different white bar.
#
# What is the main takeaway? K matters because it induces a tradeoff. It shouldn't be too large or too small (as we believe at the moment). And, generalization error estimation is very variable in low n. To see this, go back and increase n.
#
# ## Reducing variance with Cross Validation (i.e. K-fold CV)
#
# We saw previous there was a lot of variance in generalization error estimation. We can reduce some of this variance by using a very simple trick. We can rotate the train-test split so that each observation will be in the test set once. How many times is this done? K. Now we see the reason for the definition of K as it tells you how many times you validate. Why is it called "cross"? Because the training set crosses over as it does the rotation. Each observation is inside a training set K-1 times. This point will become important later. Why is it called K-fold? Because a fold is one set of training-test and there are K unique folds during the whole procedure.
#
# How does this work? Well, let's say K=10, a typical value. This means in each "fold", 90% of the data is in the training set and 10% of the data is in the test set. As we run through the K folds, we train a model on the training set and predict on the test set and compute oos residuals We aggregate those oos residuals over the folds to result in n oos residuals. We then run our error metric on all n.
#
# Let's begin with the dataset from the previous demo. Here is that code that will create the folds by specifying the K=10 test sets by index. The training sets can then be found by the set difference function.
#
# We're gonna jump right into using sci-kit learn's KFold module. 

# +
# importing dependencies
from sklearn.model_selection import KFold

# initializing model
model = LinearRegression() # intercept will fit by default

# 10 folds
K = 10
kf = KFold(n_splits=K)

# list to hold residuals
oos_cross_validation_residuals = []

for train_index , test_index in kf.split(x):
    X_train = pd.DataFrame([x[i] for i in train_index])
    y_train = pd.DataFrame([y[i] for i in train_index])

    X_test = pd.DataFrame([x[i] for i in test_index])
    y_test = pd.DataFrame([y[i] for i in test_index])
    
    model.fit(X_train,y_train)
    y_hat_g = model.predict(X_test)
     
    oos_cross_validation_residuals += [np.std(y_test - y_hat_g)]
    
oos_cross_validation_residual_average = np.std(oos_cross_validation_residuals)
oos_cross_validation_residual_average
# -

# How does this CV error look over K? 

# +
# import divisors() method from sympy
from sympy import divisors

Kuniques = list(set(divisors(n)) - set([1]))
Kuniques

# +
# constants
Nsim_per_K = 1 # simulations per K fold
num_Kuniques = len(Kuniques)
Ks = Kuniques * Nsim_per_K
n = len(y)

oos_residuals_CV = pd.DataFrame({'standard_error': [0 for i in range(Nsim_per_K * num_Kuniques)],
                        'K': [0 for i in range(Nsim_per_K * num_Kuniques)]})

for i in range(len(Ks)):
    # constants
    K = Ks[i]
    prop_train = (K - 1) / K
    n_train = int(prop_train * n)

    # train test split indices
    indices = set([i for i in range(n)])
    train_index = set(random.sample(indices, n_train))
    test_index = indices - train_index
    
    # cast as lists
    train_index = list(train_index)
    test_index = list(test_index)

    # split based on indices
    X_train = pd.DataFrame([x[i] for i in train_index])
    y_train = pd.DataFrame([y[i] for i in train_index])

    X_test = pd.DataFrame([x[i] for i in test_index])
    y_test = pd.DataFrame([y[i] for i in test_index])

    # fitting model
    g_mod = LinearRegression().fit(X_train, y_train)
    
    # get predictions
    y_hat_g = g_mod.predict(X_test)
    
    # get error metric
    if len(y_test) == 1:
        g_s_e = np.absolute([float(y_test.iloc[i]) - y_hat_g[i] for i in range(len(y_test))])
    else:
        g_s_e = np.std([float(y_test.iloc[i]) - y_hat_g[i] for i in range(len(y_test))])
    
    # update results
    oos_residuals_CV.loc[i] = np.asarray([g_s_e, K])
    
# output results
oos_residuals_CV
# -

# There is still an effect of the one random fold. Let's do this many times and look at the distribution just like before.

# +
# constants
Nsim_per_K = 500 # simulations per K fold
num_Kuniques = len(Kuniques)
Ks = Kuniques * Nsim_per_K
n = len(y)

oos_residuals = pd.DataFrame({'standard_error': [0 for i in range(Nsim_per_K * num_Kuniques)],
                        'K': [0 for i in range(Nsim_per_K * num_Kuniques)]})

for i in range(len(Ks)):
    # constants
    K = Ks[i]
    prop_train = (K - 1) / K
    n_train = int(prop_train * n)

    # train test split indices
    indices = set([i for i in range(n)])
    train_index = set(random.sample(indices, n_train))
    test_index = indices - train_index
    
    # cast as lists
    train_index = list(train_index)
    test_index = list(test_index)

    # split based on indices
    X_train = pd.DataFrame([x[i] for i in train_index])
    y_train = pd.DataFrame([y[i] for i in train_index])

    X_test = pd.DataFrame([x[i] for i in test_index])
    y_test = pd.DataFrame([y[i] for i in test_index])

    # fitting model
    g_mod = LinearRegression().fit(X_train, y_train)
    
    # get predictions
    y_hat_g = g_mod.predict(X_test)
    
    # get error metric
    if len(y_test) == 1:
        g_s_e = np.absolute([float(y_test.iloc[i]) - y_hat_g[i] for i in range(len(y_test))])
    else:
        g_s_e = np.std([float(y_test.iloc[i]) - y_hat_g[i] for i in range(len(y_test))])
    
    # update results
    oos_residuals.loc[i] = np.asarray([g_s_e, K])
    
# output results
oos_residuals
# -

# What is the variability?

# +
# mean of the standard errors grouped by K
oos_results_summary = pd.DataFrame({'K': list(oos_residuals.groupby(['K']).mean().index),
                                 'K_avg': list(oos_residuals.groupby(['K']).mean()['standard_error']),
                                 'K_SE': list(oos_residuals.groupby(['K']).std()['standard_error'])
                                }).set_index('K')

oos_results_summary

# +
sns.displot(oos_residuals,
             x = 'standard_error',
             hue = 'K',
             kind = 'kde',
             multiple="stack",
             palette = 'bright',
             alpha = 0.9,
             aspect = 1.5) # width = height * aspect

plt.axvline(label='Mean 2.0', x=oos_results_summary.iloc[0, 0], color='blue')
plt.axvline(label='Mean 5.0', x=oos_results_summary.iloc[1, 0], color='orange')
plt.axvline(label='Mean 10.0', x=oos_results_summary.iloc[2, 0], color='green')
plt.axvline(label='Mean 25.0', x=oos_results_summary.iloc[3, 0], color='red')
plt.axvline(label='Mean 50.0', x=oos_results_summary.iloc[4, 0], color='purple')
plt.axvline(label='CV', x=oos_cross_validation_residual_average, color='black')

# -

# Admittedly, I don't know the properties of CV estimates as well as I should. Thus, there will be only procedural questions on the next exam. I do know that selecting K "optimally" for general datasets is an open question.
#
# There is one other nice thing about having folds, you can estimate the standard error in your generalization estimate by pretending you have K iid samples and pretending the normal theory applies. For example, let's say K = 5. Instead of aggregating all residuals, we leave them separate and get K = 5 difference estimates for generalization error.

# +
# importing dependencies
from sklearn.model_selection import KFold

# initializing model
model = LinearRegression() # intercept will fit by default

# 10 folds
K = 10
kf = KFold(n_splits=K)

# list to hold residuals
oos_cross_validation_residuals = []

for train_index , test_index in kf.split(x):
    X_train = pd.DataFrame([x[i] for i in train_index])
    y_train = pd.DataFrame([y[i] for i in train_index])

    X_test = pd.DataFrame([x[i] for i in test_index])
    y_test = pd.DataFrame([y[i] for i in test_index])
    
    model.fit(X_train,y_train)
    y_hat_g = model.predict(X_test)
     
    oos_cross_validation_residuals += [np.std(y_test - y_hat_g)]
    
oos_cross_validation_residual_mean = np.mean(oos_cross_validation_residuals)
print(oos_cross_validation_residual_mean)

oos_cross_validation_residual_std = np.std(oos_cross_validation_residuals)
print(oos_cross_validation_residual_std)

# +
# calculating confidence interval
import scipy.stats as st

confidence_interval = st.t.interval(alpha = 0.95,
                                   df = len(y) - 1,
                                   loc = oos_cross_validation_residual_mean,
                                   scale = oos_cross_validation_residual_std)

print(confidence_interval)
print(min(y), max(y))
# -

# Although this is technically nonsense since they're not iid samples since the training set is crossed over containing mostly the same observations, at least it's something. In the above example, we've managed to capture the true generalization error.
#
# Coverage in this confidence interval is over D. So I wouldn't gain much insight by simulating different splits with the same K.
#
# I believe confidence intervals for generalization error is an open problem or maybe proved that you can't find them in general situations.
#
# Here's a real data example with the `breast_cancer` dataset.

# +
from sklearn.datasets import load_breast_cancer

data = load_breast_cancer(as_frame = True)
df = data.frame
X = df.iloc[:,:-1]
y = df.iloc[:,-1]

X

# +
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import KFold

k = 5
kf = KFold(n_splits=k, random_state=None)

model = LinearRegression()
 
residuals = []
 
for train_index , test_index in kf.split(X):
    X_train , X_test = X.iloc[train_index,:],X.iloc[test_index,:]
    y_train , y_test = y[train_index] , y[test_index]
     
    model.fit(X_train,y_train)
    y_hat_g = model.predict(X_test)
     
    residuals += [np.std(y_test - y_hat_g)]

print(np.mean(residuals))
print(np.std(residuals))

# +
# calculating confidence interval
import scipy.stats as st

confidence_interval = st.t.interval(alpha = 0.95,
                                   df = len(X) - 1,
                                   loc = np.mean(residuals),
                                   scale = np.std(residuals))

print(confidence_interval)

StarcoderdataJupyterScriptsDedupFiltered