<jupyter_start><jupyter_text>Show avg prediction and variance for synthetic data<jupyter_code>import numpy as np
import pandas as pd

from sklearn.linear_model import LinearRegression, Ridge, Lasso, LogisticRegression
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.datasets import load_boston, load_iris, load_wine, load_digits, \
                             load_breast_cancer, load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix, precision_score, recall_score

from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor
from sklearn.metrics import mean_absolute_error

import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'

from rfpimp import *

from sklearn import tree
from dtreeviz.trees import *

def jeremy_trick_RF_sample_size(n):
    # Jeremy's trick
    from sklearn.ensemble import _forest as forest
    forest._generate_sample_indices = \
        (lambda rs, n_samples: forest.check_random_state(rs).randint(0, n_samples,n))
def rent():
    # can't use too many records as model.score() will be average of many values
    # which reduces variance; we're trying to highlight not hide that variance
    df_rent = pd.read_csv("data/rent-ideal.csv")
    X = df_rent[['bedrooms','bathrooms','latitude','longitude']].values
    y = df_rent['price'].values

    # shuffle
    idx = np.random.permutation(len(X))
    X = X[idx]
    y = y[idx]
    return X, y

def testobs(X, y, nchunks = 200, ntests=5):
    n = len(X)
    chunk_size = n // nchunks

    # grab say just 5 records; too many and average r2 over test set shows too little variance
    if chunk_size<ntests:
        raise "Oops"
    X_test, y_test = X[0:ntests], y[0:ntests]
    return X_test, y_test, chunk_size
def showme(tree_sizes,all_test_errors,all_test_scores,max_features,min_samples_leaf=1,
          std_range=(0,.13),
          y_range=(0,0.3)):
    test_means = np.mean(all_test_errors, axis=0)
    test_stds = np.std(all_test_errors, axis=0)
    test_variances = np.var(all_test_errors, axis=0)
    test_r2_means = np.mean(all_test_scores, axis=0)
    test_r2_stds = np.std(all_test_scores, axis=0)
    test_r2_variances = np.var(all_test_scores, axis=0)

    fig, ax = plt.subplots(1, 1, figsize=(4,2.8))
    a = ax.plot(tree_sizes, test_r2_means, lw=.75, marker='.', label="test $R^2$")

    plt.fill_between(tree_sizes, 
                     test_r2_means - test_r2_stds,
                     test_r2_means + test_r2_stds,
                     color='gray', alpha=0.2, label="test $R^2$ $\sigma$")

    ax2 = ax.twinx()
    b = ax2.plot(tree_sizes, test_r2_stds, lw=1, label="test set stdev", c='orange')
    ax2.yaxis.label.set_color('orange')
    ax2.tick_params(axis='y', colors='orange')
    ax2.set_ylim(*std_range)

    ax.legend(loc='lower right')
    ax.set_xlabel("Num trees in forest")
    ax.set_ylabel("Test set $R^2$")
    ax.set_ylim(*y_range)
    
    if max_features==1.0:
        ax.set_title(f"BAGGING test error vs num trees, {min_samples_leaf} samples/leaf", fontsize=10)
    else:
        ax.set_title(f"RF rest error vs num trees, {min_samples_leaf} samples/leaf", fontsize=10) 
def sim_ntrees(X_train, X_test, y_train, y_test,
                max_features = 'auto',
                min_samples_leaf=1,      # lowest bias but increases variance
                ntrials = 5,
              tree_sizes=[1,2,3,4,5,6,10,15,20,30,40,50]):

    all_test_scores = np.zeros(shape=(ntrials,len(tree_sizes)))
    for i in range(ntrials):
        print(f" {i}", end='')
        test_scores = []
        for j,n_estimators in enumerate(tree_sizes):
            rf = RandomForestRegressor(n_estimators=n_estimators,
                                       min_samples_leaf=min_samples_leaf,
                                       max_features=max_features,
                                       bootstrap=True,
                                      n_jobs=-1)
            rf.fit(X_train, y_train)
            s = mean_absolute_error(y_test, rf.predict(X_test))
    #         s = rf.score(X_test, y_test)
            test_scores.append(s)

        all_test_scores[i] = test_scores
        avg = np.mean(test_scores)
        
    return all_test_scores, tree_sizes<jupyter_output><empty_output><jupyter_text>bulldozer data<jupyter_code>def bulldozer(basefeatures = ['SalesID',
                    'MachineID', 'ModelID',
                    'datasource', 'YearMade',
                    # some missing values but use anyway:
                    'auctioneerID', 'MachineHoursCurrentMeter']):
    df_bull = pd.read_feather("data/bulldozer-train.feather") # these are sorted by date

    df_bull = df_bull.fillna(0) # flip missing numeric values to zeros

    ntest = 15_000
    df_bull_test = df_bull.iloc[-ntest:] # take only last 100,000 records
    df_bull_train = df_bull.iloc[-(100_000+ntest):-ntest] # train only last 100,000 records

    X_train, y_train = df_bull_train[basefeatures], df_bull_train['SalePrice']
    X_test, y_test = df_bull_test[basefeatures], df_bull_test['SalePrice']

    return X_train, X_test, y_train, y_test<jupyter_output><empty_output><jupyter_text>OverfittingCan try playing with min_samples_leaf, but better to remove features that cause overfitting. SalesID is unique record ID never to be seen again.  Remove that and accuracy on validation set goes higher. Lesson: it's about the features not the model so much.<jupyter_code>features=[
    'MachineID', 'ModelID',
    'datasource', 'YearMade',
    'auctioneerID', 'MachineHoursCurrentMeter']

X_train, X_test, y_train, y_test = bulldozer(basefeatures=features)

#max_features = len(X_train.columns)
max_features = len(X_train.columns)
min_samples_leaf = 1
ntrials = 5

all_test_scores, tree_sizes = sim_ntrees(X_train, X_test, y_train, y_test, ntrials=ntrials,
                                         max_features=max_features,
                                         min_samples_leaf=min_samples_leaf,
                                         tree_sizes=[1,2,3,4,5,6,10,15,20,30,40,50,100])
test_means = np.mean(all_test_scores, axis=0)
test_variances = np.std(all_test_scores, axis=0)
test_variances
figsize = (4,2.8)
fig, ax = plt.subplots(1, 1, figsize=figsize)

a = ax.plot(tree_sizes, test_means, lw=.75, marker='.', label="test set $\sum |y - \hat{y}|$")

plt.fill_between(tree_sizes, 
                 test_means - test_variances*3,
                 test_means + test_variances*3,
                 color='gray', alpha=0.2)

ax.set_xlabel("Num trees in forest")
ax.set_ylabel("$\sum |y - \hat{y}|$")
ax.set_ylim(6900,10_000)

if max_features==len(X_train.columns):
    ax.set_title(f"Bulldozer bagged forest error vs num trees\n{ntrials} trials, {len(X_train)} training obs., {len(X_test)} test\n{max_features} max features, {min_samples_leaf} min_samples_leaf", fontsize=10)
else:
    ax.set_title(f"Bulldozer RF error vs num trees\n{ntrials} trials, {len(X_train)} training obs., {len(X_test)} test\n{max_features} max features, {min_samples_leaf} min_samples_leaf", fontsize=10)<jupyter_output><empty_output><jupyter_text>Independent training samples<jupyter_code>nchunks = 200
X, y = rent()
X_test, y_test, chunk_size = testobs(X, y, nchunks=nchunks)
max_features = 1.0      # use this to simulate bagging but not RF
#max_features = 'auto'
min_samples_leaf=1      # lowest bias but increases variance

tree_sizes = [1,2,3,4,5,6,10,15,20,30,40,50]

all_test_scores = np.zeros(shape=(len(tree_sizes), nchunks-1)) # first chunk is test set
for i,n_estimators in enumerate(tree_sizes):
    print(f"{n_estimators} trees: ", end='')
    test_scores = []
    for j in range(1,nchunks): # first chunk is test set; skip
#         print(j,' ', end='')
        X_train = X[j*chunk_size:j*chunk_size+chunk_size]
        y_train = y[j*chunk_size:j*chunk_size+chunk_size]
        rf = RandomForestRegressor(n_estimators=n_estimators,
                                   min_samples_leaf=min_samples_leaf,
                                   max_features=max_features,
                                   bootstrap=True)
        rf.fit(X_train, y_train)
        s = rf.score(X_test, y_test)
        test_scores.append(s)

    all_test_scores[i] = test_scores
    avg = np.mean(test_scores)
    std = np.std(test_scores)
    print(avg, std)
test_means = np.mean(all_test_scores, axis=1)
test_variances = np.std(all_test_scores, axis=1)
test_variances
figsize = (4,2.8)
fig, ax = plt.subplots(1, 1, figsize=figsize)

a = ax.plot(tree_sizes, 1-test_means, lw=.75, marker='.', label="test set 1-R^2")

plt.fill_between(tree_sizes, 
                 1-test_means - test_variances,
                 1-test_means + test_variances,
                 color='gray', alpha=0.2)

ax2 = ax.twinx()
b = ax2.plot(tree_sizes, test_variances, lw=1, label="test set $\sigma$", c='orange')
ax2.yaxis.label.set_color('orange')
#ax2.set_ylabel("Test $\sigma$")
ax2.tick_params(axis='y', colors='orange')
ax2.set_ylim(-0.002,.8)

ax.set_xlabel("Num trees in forest")
ax.set_ylabel("1-$R^2$")
ax.set_ylim(-.05,1.5)
lines = a+b
labels = [l.get_label() for l in lines]
ax.legend(lines,labels,loc='upper right')

ax.text(5,.00,f"ntrain={len(X_train)}, ntest={len(X_test)}, ntrials={nchunks}")

if max_features==1.0:
    ax.set_title("NYC Rent bagged forest error vs num trees", fontsize=10)
else:
    ax.set_title("NYC Rent RF error vs num trees", fontsize=10)<jupyter_output><empty_output><jupyter_text>Synthetic data<jupyter_code>def syndata(n = 100, ncols = 50):
    "Eqn 15.8 p599 ESLII Hastie et all"
    df = pd.DataFrame()
    for i in range(ncols):
        df[f'x{i}'] = np.random.normal(size=n)
    X = df.values
    y = np.sum(X,axis=1) / np.sqrt(ncols)# + np.random.normal(size=n)
    return X, y
def simulate(ntrials, max_features, min_samples_leaf=1, bootstrap=True):
#     tree_sizes = [1,2,3,4,5,6,7,8,9,10,20,30]#,40,50]#,200,400]
    tree_sizes = [1,2,5,8,10,15,20,30,40,50,100, 200, 300]
    all_test_errors = np.zeros(shape=(ntrials, len(tree_sizes)))
    all_test_scores = np.zeros(shape=(ntrials, len(tree_sizes)))

    # single test set
    X_test, y_test = syndata(n=100)

    print("Trial: ", end='')
    for i in range(ntrials):
        print(f"{i} ", end='')
        X_train, y_train = syndata(n=100) # new train set each trial

        error_per_tree_size = []
        r2_per_tree_size = []
        for n_estimators in tree_sizes:
            rf = RandomForestRegressor(n_estimators=n_estimators,
                                       min_samples_leaf=min_samples_leaf,
                                       max_features=max_features,
                                       bootstrap=bootstrap,
                                       n_jobs=-1)
            rf.fit(X_train, y_train)
#             y_pred = rf.predict(X_test)
#             e = mean_absolute_error(y_test, y_pred)
#             error_per_tree_size.append(e)
            s = rf.score(X_test, y_test)
            r2_per_tree_size.append(s)

#         all_test_errors[i] = error_per_tree_size
        all_test_scores[i] = r2_per_tree_size
    return all_test_errors, all_test_scores, tree_sizes
all_test_errors, all_test_scores, tree_sizes = simulate(ntrials=10, max_features='auto')
showme(tree_sizes, all_test_errors, all_test_scores, max_features='auto')
all_test_errors, all_test_scores, tree_sizes = simulate(ntrials=20, max_features=1.0)
showme(tree_sizes, all_test_errors, all_test_scores, max_features=1.0)
all_test_errors, all_test_scores, tree_sizes = simulate(ntrials=20, max_features=1)
showme(tree_sizes, all_test_errors, all_test_scores, max_features=1)
all_test_errors, all_test_scores, tree_sizes = simulate(ntrials=100, max_features='auto')
showme(tree_sizes, all_test_errors, all_test_scores, max_features='auto')
all_test_errors, all_test_scores, tree_sizes = simulate(ntrials=100, max_features=1.0)
showme(tree_sizes, all_test_errors, all_test_scores, max_features=1.0)
all_test_errors, all_test_scores, tree_sizes = simulate(ntrials=100, max_features=1)
showme(tree_sizes, all_test_errors, all_test_scores, max_features=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11691091717437839, 'fraction_numerical': 0.026181458916515726, 'mean_word_length': 3.3584499461786868, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '618687', 'license': 'MIT', 'n_tokens_mistral': 4633, 'n_tokens_neox': 4334, 'n_words': 838, 'path': 'notebooks/trees/viz-variance.ipynb', 'repo_name': 'dashbrookhart/msds621', 'size': '12146'}