<jupyter_start><jupyter_text>Learning Curve and Model Score EstimationTODO: expand to: Emperical study of model score estimation. The [King County Housing Data](https://www.kaggle.com/harlfoxem/housesalesprediction/data) from Kaggle will be used.  The model will predict house price.The Learning Curve for a particular model, is the model's score vs amount of data.  More data allows for a better trained model which allows for a higher score (or lower error).The Learning Curve can also be used to determine how much data is necessary to train on.  Once the curve flattens out, additional data is no longer necessary.In the following, RMSE of the predicted price vs the actual price will be computed.  Scikit Learn has 'neg_mean_squared_error', but RMSE is more natural and it is easy to create a custom scorer.Scikit Learn maximizes based on the scorer.  If the scorer is a measure of error, such as RMSE or MSE, then the negative of the value needs to be used so that error is minimized.The Learning Curve is created by estimating the model's score.  The most often used "estimator" of model score is K-Fold Cross Validation, but other model score estimators are possible.The model score estimator itself has bias and variance.  This can be seen in the learning curve.This notebook will create a Learning Curve and will examine the bias and variance of model score estimators.   Common Imports and Notebook Setup<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn as sk
%matplotlib inline
sns.set() # enable seaborn style<jupyter_output><empty_output><jupyter_text>Prior Notebook SetupAs per [Cross Validation: W and W/O cross_val_score](http://nbviewer.jupyter.org/github/sdiehl28/tutorial-jupyter-notebooks/blob/master/projects/king_county/CV_Pipeline.ipynb)<jupyter_code># read king county housing data
parser = lambda date: pd.datetime.strptime(date, '%Y%m%dT%H%M%S')
all_data = pd.read_csv('./data/kc_house_data.csv', parse_dates=['date'], date_parser=parser)

# compute mean sales price per house
mean_price_by_id = all_data.groupby(['id'])['price'].mean()

# convert this to a data frame
mean_price_by_id = pd.DataFrame(mean_price_by_id).reset_index()

# join on id to get a new column for mean price
all_data = all_data.merge(right=mean_price_by_id, on='id', suffixes=['_orig', '_mean'])

# drop original price and de-dup the records
all_data = all_data.drop('price_orig', axis=1).drop_duplicates(subset='id')

# remove id and date fields
all_data = all_data.drop(['id', 'date'], axis=1)

# repace yr_built with yr_renovated if yr_renovated != 0
mask = all_data['yr_renovated'] != 0
all_data.loc[mask, 'yr_built'] = all_data.loc[mask, 'yr_renovated']
all_data = all_data.drop('yr_renovated',axis=1)

all_data = all_data.drop(['zipcode'], axis=1)

X_all = all_data.drop('price_mean', axis=1)
y_all = all_data['price_mean']
print(X_all.shape)
print(y_all.shape)<jupyter_output>(21436, 16)
(21436,)<jupyter_text>Learning CurveThe full King County Housing dataset consists of 21613 records.  We can simulate having less observations, from the same distribution, by taking a random subsample of these records.  In the following the size of the sample varies from 10% to 100% of the total number of records.A given random sample of data may be unrepresentative of the data as a whole. This sample variance leads to differences in model score estimations.  To reduce the sample variance, each random sample will be chosen 10 times from the original King County dataset and the average K-Fold CV scores will be used.  Note that this is not the same thing as Repeated K-Fold CV.  Repeated K-Fold CV uses the *same* sample, and makes repeated K-Fold CV model score estimates using different train/test splits.  Repeated K-Fold CV cannot account for sample variance, but it can account, to some extent, for the inherent variance in the K-Fold CV estimation itself.With the subsample is 100% of the original data, there is only one way to chose the data.  For this case, the code below does amount to Repeated K-Fold CV.Repeated K-Fold CV has been shown to reduce the variance of the model score's estimate, but it does not cause the bias to approach zero. The modern (2018) recommendation appears to be to use 1 to 20 iterations of K-Fold CV and to consider the variance implied by repeated K-Fold CV as being too small.  That is, the true value of the model's score may be more than 2 standard deviations away from the point-estimate, more than 5% of the time, when using Repeated K-Fold CV.TODO: cite references for previous paragraph More DataAs the amount of data goes up, the RMSE and the standard deviation of RMSE will go down. Higher K in K-Fold CVFor a given amount of data, as K goes up, there is more data available for training and hence the bias of the model's score does down.  However there is less data available for validation and hence the variance of the model's score goes up. In what follows, K will be varied to observe these relationships.The flatter the learning curve is, the less affect K will have on the estimate of the model's score. Model Score Estimation is Not Model CreationThe purpose of CV is to provide an estimate of how well a specific model will do on out-of-sample data.  The estimate of the model's performance can be used for:* Model Selection* Model EstimationIn model selection, we use the estimate of the model's score to chose the best model (or the best hyperparameters). In most cases, the bais in the estimate of the model's score is similar for the models being considered.  As such, this bias does not impact the selection of which model is best.In mode estimation, we use the estimate of the model's score to estimate how well the model will do when deployed.A common misunderstanding is that CV is used to create or improve a specific model.  This is not correct.  In particular, no model created during the CV process is used for model deployment.  The model that is created for deployment is trained on the entire dataset (unless the learning curve shows that using all the data is unnecessary).GridSeachCV() and RandomSearchCV() appear to use CV to create the best model, but actually there is a default parameter of "refit=True", which invokes a final step of refitting the model using all available data to the best hyperparameters found. In this case too, no model created during CV is returned to the user.As can be seen from the Learning Curve, CV provides a pessimistic estimate of the model's performance.  CV uses less than the full amount of data to train with, and if the learning curve is steep, then this impacts the CV estimate of the model's score significantly.  If the Learning Curve is flat, then CV will not have a pessimistic bias of the model's performance.It is not always best to pick the model with the best score.  If the estimate of the best score has high variance, then the selected model might not really be best.  ESL recommends the one sigma rule, which is to pick the simplest model that is within one standard deviation of the best model.  Ideas such as this will be explored in greater detail in the notebook on Model Selection and Bias Correction.TODO: Link to ESL (and ISL if the above is stated there).All of the above is graphically demonstrated below.   Imports<jupyter_code>from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.neighbors import KNeighborsRegressor
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, KFold, RepeatedKFold
from sklearn.metrics import mean_squared_error
from time import time<jupyter_output><empty_output><jupyter_text>DaskDask offers at least two advantages:* ability to work with problems that do not fit in a single computer's memory* distributed processing across computers (and excellent use of all cores on a single computer)Below, Dask is use for distributed processing in order to speed up the computations.Scikit Learn has built-in multiprocessing support for at least the following:* cross_val_score()* GridSearchCV()* RandomSearchCV()This can be enabled by setting n_jobs=-1, for example.In my performance tests, using Dask improved performance even on a single computer.  Using Dask for distributed computing across two computers offered a major performance improvement.If you are using more than one computer, you should ensure that your Anaconda environment is identical on all computers.My setup consists of two desktop computers:* Intel Xeon E5-1630v3 3.70GHz -- 4 physical cores with hyperthreading* Intel Core I7-7700 3.60GHz -- 4 physical cores with hyperthreading<jupyter_code># Dask Setup
# import joblib

from dask.distributed import Client
try:
    client.close()
except:
    print('client was not open')
client = Client('atlantis:8786')
client<jupyter_output>client was not open<jupyter_text>Compute Learning Curve10 different realizations each of 10 different sized subsets of the data will be used.  K will be varied to see the effect of K on model score estimation.<jupyter_code># Create custom scorer for RMSE
from sklearn.metrics import make_scorer

def my_rmse_scorer(y, y_pred):
    
    # use numpy vectorized arthimetic
    return -np.sqrt(np.mean((y-y_pred)*(y-y_pred)))

rmse_scorer = make_scorer(my_rmse_scorer)
# Create random subsample
def get_random_subsample(X, y, frac, random_state=None):
    ''' 
    X is a 2D numpy array of numeric features
    y is a 1D numpy array of numeric target values
    frac is the percentage of the sample to take
    returns a subsample of X, y as a tuple
    '''
    
    # affects numpy.random.choice
    if random_state:
        np.random.seed(random_state)
    
    # select a random sample of frac size
    n_records = X.shape[0]
    n = int(frac*n_records)
    idx = np.random.choice(n_records, size=n, replace=False)
    X = X[idx, :]
    y = y[idx]
    
    return (X, y)
# define a specific model
def create_model():
    # The hyperparameters used here were determined in the CV_Hyperparameter notebook
    scale = StandardScaler()
    pca = PCA(n_components=8)
    knn = KNeighborsRegressor(n_neighbors=10, weights='distance')

    return Pipeline(steps=[('scale', scale),('pca',pca), ('knn',knn)])
# Learning Curve4
def learning_curve4(classifier, X_s, y_s, k_fold, n_repeats, random_state=None):
        
    # print(f'RandomState: {random_state} frac: {frac:3.1f} k_fold {k_fold}')

    cv = RepeatedKFold(n_splits=k_fold, n_repeats=n_repeats, random_state=random_state+9)

    # compute the K cross validated scores
    # use n_jobs = 1 as this entire function will be parallelized
    scores = cross_val_score(classifier, X_s, y_s, cv=cv, n_jobs=1, scoring = rmse_scorer)

    return pd.DataFrame(data=[[k_fold, X_s.shape[0], random_state, scores.mean(), scores.std()]], 
                            columns=['k_fold', 'len', 'random_state', 'cv_score_mean', 'cv_score_std'])
# Dask Performance Optimization for my setup only

# Dask scheduler distributes jobs to computers based on number of cores
# Both my computers have same number of cores, but the I7 computer (.67) is faster
w = [worker for worker in client.scheduler_info()['workers'].keys()]
if '67' in w[0]:
    w_fast = w[0]
    w_slow = w[1]
else:
    w_fast = w[1]
    w_slow = w[0]
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

# Compute Learning Curve Dask submit
start = time()

classifier = create_model()

r = []
repeated_cv = 1
for random_state in range(20):
    for frac in np.linspace(0.1, 1.0, 10):
        X_s, y_s = get_random_subsample(X_all.values, y_all.values, frac, random_state)
        for k_fold in [2, 3, 5, 10, 20]:
            r.append(client.submit(learning_curve4, classifier, X_s, y_s, k_fold, repeated_cv, random_state))

# gather the results from the scheduler
dfs_20_1 = client.gather(r)

print(f'Total Time: {time() - start}')
df = pd.concat(dfs_20_1, ignore_index=True)
df.head()
# group by K and Number of Records in Preparation
# for plotting the Learning Curve
rmse_data = df.groupby(['k_fold', 'len'])['cv_score_mean'].size()
rmse_data.head(11)<jupyter_output><empty_output><jupyter_text>We see there are 20 values, cooresponding to the 20 realizations of a subset of the full data, for which K-Fold Cross Validation was performed.  The average of these values will be used for the Learning Curve.<jupyter_code>rmse_data = df.groupby(['k_fold', 'len'])['cv_score_mean'].mean()
rmse_data.head(11)
# plot RMSE vs Sample Size for various values of K
plt.figure(figsize=(12,6))

rmse_data = df.groupby(['k_fold','len'])['cv_score_mean'].mean()
for k in [2, 3, 5, 10, 20]:
    plt.plot(rmse_data.loc[k].index.values, -rmse_data.loc[k].values, marker='o', label=k)

plt.legend(title="K Folds")
plt.xlabel('Records')
plt.ylabel('RMSE')
plt.ylim(150000, 210000)
plt.title('Learning Curve -- 20 Samples')<jupyter_output><empty_output><jupyter_text>Learning Curve AnalysisPerhaps the first thing to note is that the general shape of the above Learning Curve is downward sloping and the shape is about the same regardless of the value of K.For determining how much data is necessary, wherever the curve levels off, regardless of K, is good.  The Learning Curve appears to be leveling off at about the 90% mark, which is 19292 records.K=20 provided the lowest RMSE score at every data set size.  Does this mean that we should choose K=20 to create a better model? The short answer is "no".The RMSE scores that are plotted are estimates of the model's actual score.  These estimates have no affect on the model's true score.The reason K=20 produces a lower estimated RMSE is that is has less pessimistic bias.  It has less pessimistic bias because it is using more of the data and the Learning Curve is not flat.The difference between each learning curve, for a given dataset size, is the difference in bias due to having less data to train on. Sample VarianceThe sample variance in the above plot was mitigated by using 10 different realizations of each subset of the full data.  This hides the variance of higher K values.  Below the same Learning Curve will be plotted using just 1 realization of each subset of data.************Because the bias drops with higher K, leave-one-out CV (the highest K possible) is sometimes used to get the least biased estimate of the model's performance.  However although it is least biased, it has the most variance, as seen below.If K-Fold CV is being used for Model Selection, and the lowest CV score is used, the wrong model may be selected if the variance is too high!<jupyter_code># insted of taking the mean of the 10 values, just take the first value
dfg = df.groupby(['k_fold','len'])['cv_score_mean']
rmse_1 = dfg.apply(lambda x: x.iloc[0])
rmse_1.head(11)
# plot RMSE vs Sample Size for various values of K
plt.figure(figsize=(12,6))

for k in [2, 3, 5, 10, 20]:
    plt.plot(rmse_1.loc[k].index.values, -rmse_1.loc[k].values, marker='o', label=k)

plt.legend(title="K Folds")
plt.xlabel('Records')
plt.ylabel('RMSE')
plt.ylim(150000, 210000)
plt.title('Learning Curve -- Single Subsample CV')<jupyter_output><empty_output><jupyter_text>The variance in the estimate of the model accuracy is too high to produce a meaningful Learning Curve.  This is what Repeated K-Fold CV was designed to help alleviate.Run Repeated K-Fold CV on a one sample of data at each size.<jupyter_code>import warnings
warnings.filterwarnings("ignore", category=UserWarning)

# Compute Learning Curve Dask submit
start = time()

classifier = create_model()

r = []
repeated_cv = 20
for random_state in range(1):
    for frac in np.linspace(0.1, 1.0, 10):
        X_s, y_s = get_random_subsample(X_all.values, y_all.values, frac, random_state)
        for k_fold in [2, 3, 5, 10, 20]:
            r.append(client.submit(learning_curve4, classifier, X_s, y_s, k_fold, repeated_cv, random_state))

# gather the results from the scheduler
dfs_1_20 = client.gather(r)

print(f'Total Time: {time() - start}')
df = pd.concat(dfs_1_20, ignore_index=True)
df.head()
# plot RMSE vs Sample Size for various values of K
plt.figure(figsize=(12,6))

rmse_data = df.groupby(['k_fold','len'])['cv_score_mean'].mean()
for k in [2, 3, 5, 10, 20]:
    plt.plot(rmse_data.loc[k].index.values, -rmse_data.loc[k].values, marker='o', label=k)

plt.legend(title="K Folds")
plt.xlabel('Records')
plt.ylabel('RMSE')
plt.ylim(150000, 210000)
plt.title('Learning Curve -- Single Sample 20 Repeated CV')<jupyter_output><empty_output><jupyter_text>Repeated K-Fold CV does not make up for sample variance (the variance due to which random sample was chosen).  There is less sample variance at higher percentages of the data.At about 60% of the data (12681 records) and above, the Learning Curve looks reasonable.Normally sample variance cannot be avoided as we only have one sample of the data.<jupyter_code>def my_cv(B, size, random_state=None):
    if random_state:
        np.random.seed(random_state)
        
    cv = []
    for i in range(B):
        idx = np.random.choice(size, size=size, replace=True)
        idx_all = set(range(size))
        idx_validate = idx_all - set(idx)
        idx_validate = np.array(list(idx_validate))
        cv.append((idx, idx_validate))
    return cv
# Learning Curve5
def learning_curve5(classifier, X, y, B, random_state=None):
        
    # print(f'RandomState: {random_state} frac: {frac:3.1f} k_fold {k_fold}')
    
    # B bootstrap samples
    cv = my_cv(B, len(X_s), random_state)

    # compute the K cross validated scores
    # use n_jobs = 1 as this entire function will be parallelized
    scores = cross_val_score(classifier, X_s, y_s, cv=cv, n_jobs=1, scoring = rmse_scorer)

    return pd.DataFrame(data=[[X_s.shape[0], random_state, scores.mean(), 
                               scores.std(), B, scores]], 
                            columns=['len', 'random_state', 'cv_score_mean', 
                                     'cv_score_std', 'B', 'cv_scores'])
# Compute Learning Curve Dask ML
B = 100
classifier = create_model()
start = time()

# let Dask ML take care of cross_val_score()
dfs_bs = []
for random_state in range(6):
    for frac in np.linspace(0.1, 1.0, 10):
        X_s, y_s = get_random_subsample(X_all.values, y_all.values, frac, random_state)
        dfs_bs.append(client.submit(learning_curve5, classifier, X_s, y_s, B, random_state))
    
dfs_bs = client.gather(dfs_bs)    

print(f'Total Time: {time() - start}')
dfs_bs = client.gather(dfs_bs)
df_bs = pd.concat(dfs_bs, ignore_index=True)
df_bs.head()
# messy code to get a numpy array from an element in a dataframe
from collections import defaultdict
scores = defaultdict(list)
for key, group in df_bs.groupby(['len','random_state'])['cv_scores']:
    scores[key[0]].extend(group.values[0])
x=[]
y=[]
for key in scores.keys():
    x.append(key)
    y.append(-np.array(scores[key]).mean())
    print(f'{key: 6d} {-np.array(scores[key]).mean(): 9.1f}')
y
# Compute Learning Curve Dask ML

random_state = 1
B = 600
classifier = create_model()
start = time()

# let Dask ML take care of cross_val_score()
dfs_bs = []
with parallel_backend('dask.distributed', scheduler_host='atlantis:8786'):
    for frac in np.linspace(0.1, 1.0, 10):
        X_s, y_s = get_random_subsample(X_all.values, y_all.values, frac, random_state)
        dfs_bs.append(learning_curve5(classifier, X_s, y_s, B, random_state))

print(f'Total Time: {time() - start}')
plt.plot(x,y)
# plot RMSE vs Sample Size for various values of K
plt.figure(figsize=(12,6))

plt.plot(x,y, marker='o', label='600-bs')
# plt.plot(rmse_data.loc[2].index.values, -rmse_data.loc[2].values, marker='o', label='20CV2')

# plt.legend(title="Bootstrap")
plt.xlabel('Records')
plt.ylabel('RMSE')
plt.ylim(150000, 230000)
# plt.legend()
plt.title('Learning Curve -- Bootstrap 600 Replicas 1 Sample')
df_bs = pd.concat(dfs_bs, ignore_index=True)
# plot RMSE vs Sample Size for various values of K
plt.figure(figsize=(12,6))

plt.plot(df_bs['len'], -df_bs['cv_score_mean'], marker='o', label='600-bs')
plt.plot(rmse_data.loc[2].index.values, -rmse_data.loc[2].values, marker='o', label='20CV2')

# plt.legend(title="Bootstrap")
plt.xlabel('Records')
plt.ylabel('RMSE')
plt.ylim(150000, 210000)
plt.legend()
plt.title('Learning Curve -- Bootstrap 600 Replicas 1 Sample')<jupyter_output><empty_output><jupyter_text>As more data becomes available, the model has more data to train on so the model is more stable, and the model has more data to be validate against, so the validation is more stable, hence the standard deviation of the RMSE goes down.For a given dataset size, as K becomes lower, more data is available for validation.  With more data available for validation, the sample variance of the validation set goes down, allowing for a more stable estimate of the model's score.  Which Value of K is "Best"The Bias-Variance tradeoff means that for a more accurate estimate of the model's performance, the higher the variance of that estimate.  As such, there is no "best" value.For the above model, it might be argued that K=5 results in a good tradeoff of bias vs variance.  Bootstrap EstimateThe bootstrap estimate of the model's score is not often used as the result is similar to 2-Fold CV but requires more computation.  The bootstrap replica for training contains about 62.8% unique records and the bootstrap replica for testing is the 37.2% of the records not in training.Scikit Learn does not have a bootstrap cross validation generator, but it is easy to create one.<jupyter_code>from scipy import stats
plt.figure(figsize=(6,4), dpi=120)

# stats.boxcox()
             
for key in scores.keys():
    if key > 10000:
        sns.distplot(-np.array(scores[key]),hist=False, label=key)
plt.legend()
t = -np.array(scores[12861])
%%R
library(boot)

bs_mean <- function(data, indices) {
  d <- data[indices] # allows boot to select sample 
  return(mean(d))
}

set.seed(108)
r_boot_ci <- function(x, t, bs_all_stat) {
  results <- NULL
  results$data <- x
  results$t <- t
  results$t0 <- bs_all_stat
  
  results$R <- length(t)  
  results$strata <- rep(1,length(x))
  results$weights <- rep(1.0/length(x), length(x))
    
  results$sim <- 'ordinary'
  results$statistic <- NULL
  results$call <- NULL
  results$stype <- 'i'
  results$seed <- get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
  
  boot.ci(results)
}
%load_ext rpy2.ipython

import pandas as pd
df = pd.DataFrame({
    'cups_of_coffee': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
    'productivity': [2, 5, 6, 8, 9, 8, 0, 1, 0, -1]
})
%%R -i df -w 4 -h 2 --units in -r 200
library(ggplot2)
ggplot(df, aes(x=cups_of_coffee, y=productivity)) + geom_line()
X = X_all.values
y = y_all.values
%%R -i X -i y
library(MASS)
print(dim(X))
print(y[20])
print(typeof(y))
boxcox(y)
from scipy.stats import boxcox
boxcox(y_all.values)
y_all.values.shape[0]
y_mean = y_all.mean()
y = y_all
N = len(y_all)
np.log(np.sum((y - y_mean)**2. / N, axis=0))
y_all.describe()/1e4
y_all.max()
y_all.mean()
plt.figure(figsize=(12,6))
sns.distplot(np.log(y))<jupyter_output>/home/agni/anaconda3/envs/anaconda/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "
{'fraction_non_alphanumeric': 0.09206294755799493, 'fraction_numerical': 0.017966639509455, 'mean_word_length': 4.488820899035067, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 39, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '335383', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7699, 'n_tokens_neox': 7180, 'n_words': 2884, 'path': 'projects/king_county_house/LC-Submit.ipynb', 'repo_name': 'sdiehl28/tutorial-jupyter-notebooks', 'size': '23321'}