<jupyter_start><jupyter_text>Tuning HyperparametersCiaran LowellMidterm ProjectUnfortunately Professor Vasiliu, I wasn't able to debug this code in time. I'm just submitting what I had so far so I can get a low partial grade instead of a 0.This project will demonstrate hyperparamaters and creation of polynomial features for the following regressors:* Lasso* Ridge* Elastic Net* SCAD* Square Root Lasso * Stepwise* Kernel* Random Forest* XGBoost<jupyter_code>import pandas as pd
#import numpy as np
import matplotlib.pyplot as plt
import keras
import xgboost as xgb
from sklearn.ensemble import RandomForestRegressor
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from sklearn.metrics import r2_score
from yellowbrick.regressor import ResidualsPlot
import sklearn.datasets
import sklearn
import sklearn.datasets as dt
from sklearn.model_selection import GridSearchCV, cross_val_score 
from sklearn.linear_model import Lasso, LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import KFold
import yellowbrick as yb
#!pip install numpy==1.15.4
import numpy as np
def DoKFold(X,y,alpha):
  PE = []
  for idxtrain, idxtest in kf.split(X):
    X_train = X[idxtrain,:]
    y_train = y[idxtrain]
    X_test  = X[idxtest,:]
    y_test  = y[idxtest]
    model = sm.OLS(y_train,X_train)
    result = model.fit_regularized(method='sqrt_lasso', alpha=alpha)
    yhat_test = result.predict(X_test)
    PE.append(MAE(y_test,yhat_test))
  return np.mean(PE)
def cv_comparison(models, X, y, cv, model_names):
    # Initiate a DataFrame for the averages and a list for all measures
    cv_accuracies = pd.DataFrame()
    count = 0
    maes = []
    mses = []
    r2s = []
    accs = []
    # Loop through the models, run a CV, add the average scores to the DataFrame and the scores of 
    # all CVs to the list
    for model in models:
        mae = -np.round(cross_val_score(model, X, y, scoring='neg_mean_absolute_error', cv=cv), 4)
        maes.append(mae)
        mae_avg = round(mae.mean(), 4)
        mse = -np.round(cross_val_score(model, X, y, scoring='neg_mean_squared_error', cv=cv), 4)
        mses.append(mse)
        mse_avg = round(mse.mean(), 4)
        r2 = np.round(cross_val_score(model, X, y, scoring='r2', cv=cv), 4)
        r2s.append(r2)
        r2_avg = round(r2.mean(), 4)
        acc = np.round((100 - (100 * (mae * len(X))) / y.sum()[0]), 4)
        accs.append(acc)
        acc_avg = round(acc.mean(), 4)
        cv_accuracies[model_names[count]] = [mae_avg, mse_avg, r2_avg, acc_avg]
        count = count + 1
    cv_accuracies.index = ['Mean Absolute Error', 'Mean Squared Error', 'R^2', 'Accuracy']
    return cv_accuracies, maes, mses, r2s, accs
from google.colab import drive
drive.mount('/content/drive')
cd /content/drive/MyDrive/Colab Notebooks/VasiliuLectures<jupyter_output><empty_output><jupyter_text>DatasetThis project will be analyzing the [Concrete Compressive Strength dataset](http://archive.ics.uci.edu/ml/datasets/concrete+compressive+strength). This dataset provides 7 input variables (amounts of elements of different ingredients in the concrete, as well as Age) that contribute to the output variable, the concrete's compressive strength.The dataset has 1030 rows.<jupyter_code>concreteDF = pd.read_csv('Concrete_DataSheet1.csv', header = 0, names = ['Cement', 'Blast Furnace Slag', 'Fly Ash', 'Water', 'Superplasticizer', 'Coarse Aggregate', 'Fine Aggregate', 'Age', 'Concrete Compressive Strength'])
concreteDF
concreteX = concreteDF.loc[:, concreteDF.columns != 'Concrete Compressive Strength']
concreteY = concreteDF.loc[:, concreteDF.columns == 'Concrete Compressive Strength']
concreteY<jupyter_output><empty_output><jupyter_text>Data PreparationBefore we begin analyzing the data, we'll scale it by the z-scores. This is crucial for some of the regressors that we'll be implementing.Rather than using a Train/Test split, we'll be evaluating all of our regressors via K-Fold analysis.<jupyter_code>scaler = sklearn.preprocessing.StandardScaler()
ConcreteX = scaler(ConcreteX)<jupyter_output><empty_output><jupyter_text>Creating Polynomial FeaturesThe Concrete Compressive Strength dataset is a complex dataset, but several of the regressors we are fitting are linear regressions.To improve their performance, before fitting them to the data, we will create polynomial features from the input data. I arbitrarily chose to make these features of the 3rd order.<jupyter_code>from sklearn.preprocessing import PolynomialFeatures

poly = PolynomialFeatures(3)
concreteXpoly = poly.fit_transform(concreteX)
concreteXpoly.shape
concreteY.dtypes<jupyter_output><empty_output><jupyter_text>LASSO Regression<jupyter_code>lasso_alpha = [int(x) for x in np.linspace(1, 100, 100)]
lasso_params = {'alpha':lasso_alpha}
lasso = GridSearchCV(sklearn.linear_model.Lasso(), param_grid = lasso_params, cv = 5, verbose = 1)
lasso.fit(concreteXpoly, concreteY)
lassotuned = lasso.best_estimator_
lasso.best_params_<jupyter_output><empty_output><jupyter_text>Ridge RegressionOur hyperparameter *alpha* is tuned to be 100,000. This value will be stored for later comparison between the regressors.<jupyter_code>ridge_alpha = [float(x) for x in np.linspace(30000, 100000, 1000)]
print(ridge_alpha)
ridge_params = {'alpha':ridge_alpha}
ridge = GridSearchCV(sklearn.linear_model.Ridge(), param_grid = ridge_params, cv = 5, verbose = 1)
ridge.fit(concreteXpoly, concreteY)
ridgetuned = ridge.best_estimator_
ridge.best_params_<jupyter_output>[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.<jupyter_text>Elastic Net<jupyter_code>eNet_alpha = [float(x) for x in np.linspace(.1, 100, 20)]
eNet_l1 = [float(x) for x in np.linspace(0.01, 1, 100)]
eNet_params = {'alpha':eNet_alpha, 'l1_ratio':eNet_l1}
eNet = GridSearchCV(sklearn.linear_model.ElasticNet(), param_grid = eNet_params, cv = 5, verbose = 0)
eNet.fit(concreteXpoly, concreteY)
eNettuned = eNet.best_estimator_
eNet.best_params_<jupyter_output>[1;30;43mStreaming output truncated to the last 5000 lines.[0m
/usr/local/lib/python3.7/dist-packages/sklearn/linear_model/_coordinate_descent.py:476: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 11011.1571758627, tolerance: 23.683356112220874
  positive)
/usr/local/lib/python3.7/dist-packages/sklearn/linear_model/_coordinate_descent.py:476: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 12625.16220039179, tolerance: 23.73902658979369
  positive)
/usr/local/lib/python3.7/dist-packages/sklearn/linear_model/_coordinate_descent.py:476: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Duality gap: 13749.371375192932, tolerance: 22.401226859563106
  positive)
/usr/local/lib/python3.7/dist-packages/sklearn/linear_model/_coordinate_descent.py:476: ConvergenceWarning: Objective did not converge. You might w[...]<jupyter_text>SCAD<jupyter_code>def scad_penalty(beta_hat, lambda_val, a_val):
    is_linear = (np.abs(beta_hat) <= lambda_val)
    is_quadratic = np.logical_and(lambda_val < np.abs(beta_hat), np.abs(beta_hat) <= a_val * lambda_val)
    is_constant = (a_val * lambda_val) < np.abs(beta_hat)
    
    linear_part = lambda_val * np.abs(beta_hat) * is_linear
    quadratic_part = (2 * a_val * lambda_val * np.abs(beta_hat) - beta_hat**2 - lambda_val**2) / (2 * (a_val - 1)) * is_quadratic
    constant_part = (lambda_val**2 * (a_val + 1)) / 2 * is_constant
    return linear_part + quadratic_part + constant_part
    
def scad_derivative(beta_hat, lambda_val, a_val):
    return lambda_val * ((beta_hat <= lambda_val) + (a_val * lambda_val - beta_hat)*((a_val * lambda_val - beta_hat) > 0) / ((a_val - 1) * lambda_val) * (beta_hat > lambda_val))
def scad_model(X,y,lam,a):
  n = X.shape[0]
  p = X.shape[1]
  # we add an extra column of 1 for the intercept
  X = np.c_[np.ones((n,1)),X]
  def scad(beta):
    beta = beta.flatten()
    beta = beta.reshape(-1,1)
    n = len(y)
    return 1/n*np.sum((y-X.dot(beta))**2) + np.sum(scad_penalty(beta,lam,a))
  
  def dscad(beta):
    beta = beta.flatten()
    beta = beta.reshape(-1,1)
    n = len(y)
    return np.array(-2/n*np.transpose(X).dot(y-X.dot(beta))+scad_derivative(beta,lam,a)).flatten()
  b0 = np.ones((p+1,1))
  output = minimize(scad, b0, method='L-BFGS-B', jac=dscad,options={'gtol': 1e-8, 'maxiter': 1e7,'maxls': 25,'disp': True})
  return output.x
def scad_predict(X,y,lam,a):
  beta_scad = scad_model(X,y,lam,a)
  n = X.shape[0]
  p = X.shape[1]
  # we add an extra column of 1 for the intercept
  X = np.c_[np.ones((n,1)),X]
  return X.dot(beta_scad)
def DoKFoldScad(X,y,lam,a):
  PE = []
  for idxtrain, idxtest in kf.split(X):
    X_train = X.iloc[idxtrain,:]
    y_train = y.iloc[idxtrain]
    X_test  = X.iloc[idxtest,:]
    y_test  = y.iloc[idxtest]
    beta_scad = scad_model(X_train,y_train,lam,a)
    n = X_test.shape[0]
    p = X_test.shape[1]
    # we add an extra columns of 1 for the intercept
    X1_test = np.c_[np.ones((n,1)),X_test]
    yhat_scad = X1_test.dot(beta_scad)
    PE.append(MAE(y_test,yhat_scad))
  return 1000*np.mean(PE)

def ScadCV(hparam):
  PEV = []
  for i in prange(len(hparam[:,0])):
    for j in prange(len(hparam[:,1])):
      PE  = []
      for idxtrain, idxtest in kf.split(X):
        X_train = X.iloc[idxtrain,:]
        y_train = y.iloc[idxtrain]
        X_test  = X.iloc[idxtest,:]
        y_test  = y.iloc[idxtest]
        beta_scad = scad_model(X_train,y_train,hparam[i,0],hparam[j,1])
        n = X_test.shape[0]
        p = X_test.shape[1]
        # we add an extra columns of 1 for the intercept
        X1_test = np.c_[np.ones((n,1)),X_test]
        yhat_scad = X1_test.dot(beta_scad)
        PE.append(MAE(y_test,yhat_scad))
    PEV.append(np.mean(PE))
  return np.array(PEV)
from numba import jit, prange
import scipy
from scipy.optimize import minimize
lam = 0.5
a = 1.001
p = concreteXpoly.shape[1]
b0 = np.random.normal(1,1,p)
y = concreteY
X = concreteX
#y = y.reshape((y.shape[0]), 1)
DoKFoldScad(X,y,0.1,1.2)<jupyter_output><empty_output><jupyter_text>Square Root Lasso<jupyter_code>#! pip install --upgrade Cython
#! pip install --upgrade git+https://github.com/statsmodels/statsmodels
import statsmodels.api as sm
from numba import jit, prange
from sklearn.metrics import mean_absolute_error as MAE
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import KFold
kf = KFold(n_splits=5,shuffle=True,random_state=1234)
def DoKFoldSqrtLasso(X,y,alpha):
  PE = []
  for idxtrain, idxtest in kf.split(X):
    X_train = X.iloc[idxtrain,:]
    y_train = y.iloc[idxtrain]
    X_test  = X.iloc[idxtest,:]
    y_test  = y.iloc[idxtest]
    model = sm.OLS(y_train,X_train)
    result = model.fit_regularized(method='sqrt_lasso', alpha=alpha)
    yhat_test = result.predict(X_test)
    PE.append(MAE(y_test,yhat_test))
  return 1000*np.mean(PE)

def SqrtLassoCV(alpha):
  PEV = []
  print(X.shape)
  print(y.shape)
  for i in prange(len(alpha[:,0])):
    print("first loop")
    PE  = []
    a = alpha[i,0]
    print(a)
    for idxtrain, idxtest in kf.split(X):
      print("inner loop")
      print(idxtrain.shape)
      print(idxtest.shape)
      X_train = X[idxtrain,:]
      print("X_train.shape = ",X_train.shape)
      y_train = y[idxtrain]
      print("y_train.shape = ", y_train.shape)
      X_test  = X[idxtest,:]
      print("X_test.shape = ", X_test.shape)
      y_test  = y[idxtest]
      model = sm.OLS(y_train,X_train)
      result = model.fit_regularized(method='sqrt_lasso', alpha=a)
      yhat_test = result.predict(X_test)
      PE.append(MAE(y_test,yhat_test))
    PEV.append(np.mean(PE))
  return np.array(PEV)
def sqrtLasso(X,y,alpha, l1_wt):
    model = sm.OLS(y,X)
    result = model.fit_regularized(method='sqrt_lasso', alpha=alpha, l1_wt = l1_wt,profile_scale=True)
    return result.predict(X), result.params
! pip install pyswarms
import pyswarms as ps
alphas = np.arange(0,10,0.1)

x_max = 20 *np.ones(1)
x_min = 0.01 * np.ones(1)
bounds = (x_min, x_max)
options = {'c1': 0.5, 'c2': 0.3, 'w': 0.9}

optimizer = ps.single.GlobalBestPSO(n_particles=20, dimensions=1, options=options, bounds=bounds)

cost, pos = optimizer.optimize(SqrtLassoCV, iters=50)<jupyter_output>2021-05-18 02:20:47,587 - pyswarms.single.global_best - INFO - Optimize for 50 iters with {'c1': 0.5, 'c2': 0.3, 'w': 0.9}

pyswarms.single.global_best:   0%|          |0/50[A<jupyter_text>Stepwise Regression<jupyter_code># Implementation of stepwise regression
def stepwise_selection(X, y, 
                       initial_list=[], 
                       threshold_in=0.01, 
                       threshold_out = 0.05, 
                       verbose=True):
    
    included = list(initial_list)
    while True:
        changed=False
        # forward step
        excluded = list(set(X.columns)-set(included))
        new_pval = pd.Series(index=excluded)
        for new_column in excluded:
            model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included+[new_column]]))).fit()
            new_pval[new_column] = model.pvalues[new_column]
        best_pval = new_pval.min()
        if best_pval < threshold_in:
            best_feature = new_pval.idxmin()
            included.append(best_feature)
            changed=True
            if verbose:
                print('Add  {:30} with p-value {:.6}'.format(best_feature, best_pval))

        # backward step
        model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()
        # use all coefs except intercept
        pvalues = model.pvalues.iloc[1:]
        worst_pval = pvalues.max() # null if pvalues is empty
        if worst_pval > threshold_out:
            changed=True
            worst_feature = pvalues.idxmax()
            included.remove(worst_feature)
            if verbose:
                print('Drop {:30} with p-value {:.6}'.format(worst_feature, worst_pval))
        if not changed:
            break
    return included<jupyter_output><empty_output><jupyter_text>Kernel Regression<jupyter_code><jupyter_output><empty_output><jupyter_text>Random Forest<jupyter_code># Number of trees in Random Forest
rf_n_estimators = [int(x) for x in np.linspace(200, 1000, 5)]
rf_n_estimators.append(1500)

# Maximum number of levels in tree
rf_max_depth = [int(x) for x in np.linspace(5, 55, 11)]
# Add the default as a possible value
rf_max_depth.append(None)

# Number of features to consider at every split
rf_max_features = ['auto', 'sqrt', 'log2']

# Criterion to split on
rf_criterion = ['mse', 'mae']

# Minimum number of samples required to split a node
rf_min_samples_split = [int(x) for x in np.linspace(2, 10, 9)]

# Minimum decrease in impurity required for split to happen
rf_min_impurity_decrease = [0.0, 0.05, 0.1]

# Method of selecting samples for training each tree
rf_bootstrap = [True, False]

# Create the grid
rf_grid = {'n_estimators': rf_n_estimators,
               'max_depth': rf_max_depth,
               'max_features': rf_max_features,
               'criterion': rf_criterion,
               'min_samples_split': rf_min_samples_split,
               'min_impurity_decrease': rf_min_impurity_decrease,
               'bootstrap': rf_bootstrap}
from sklearn.model_selection import RandomizedSearchCV
rf_base = RandomForestRegressor()

# Create the random search Random Forest
rf_tune = RandomizedSearchCV(estimator = rf_base, param_distributions = rf_grid, n_iter = 50, cv = 5, verbose = 2,  
                               n_jobs = -1, random_state = 1693)

# Fit the random search model
rf_tune.fit(concreteX, concreteY)

# View the best parameters from the random search
rf_tune.best_params_
rf_tuned = rf_tune.best_estimator_<jupyter_output><empty_output><jupyter_text>XGBoost<jupyter_code># eta (learning rate)
xgb_lr = [x for x in np.linspace(.01, .5, 100)]

#min_child_weight
xgb_mcw = [int(x) for x in np.linspace(1, 5, 5)]

# gamma (minimum loss reduction for split)
xgb_gamma = [int(x) for x in np.linspace(1,100, 10)]

# Create the grid
xgb_grid = {'eta': xgb_lr,
               'min_child_weight': xgb_mcw,
               'gamma': xgb_gamma}
xgb_base = xgb_reg
# Create the random search Random Forest
xgb_tuner = RandomizedSearchCV(estimator = xgb_base, param_distributions = xgb_grid, n_iter = 100, cv = 5, verbose = 1,  
                               n_jobs = -1, random_state = 1693)
xgb_tuner.fit(concreteX, concreteY)


# View the best parameters from the random search
xgb_tuned = xgb_tuner.best_estimator_
xgb_tuner.best_params_<jupyter_output><empty_output><jupyter_text>Neural Network<jupyter_code>from keras.models import Sequential
from keras.optimizers import Adam
from keras.callbacks import EarlyStopping
from tensorflow import keras
modelnn = Sequential()
modelnn.add(Dense(128, activation="relu"))
modelnn.add(Dense(64, activation = "relu"))
modelnn.add(Dense(32, activation="relu"))
modelnn.add(Dense(8, activation="relu"))
modelnn.add(Dense(1, activation="linear"))

modelnn.compile(loss='mean_squared_error', optimizer=Adam(lr=1e-3, decay=1e-3 / 200), metrics = [keras.metrics.MeanAbsoluteError(), keras.metrics.MeanSquaredError(), keras.metrics.MeanAbsolutePercentageError()])

es = EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=1000)

history = modelnn.fit(concreteX, concreteY, validation_split=0.3, epochs=1000, batch_size=100, verbose=0, callbacks=[es])


modelnn.evaluate(concreteX, concreteY, verbose = 3)
len(history.history['mean_squared_error'])
np.mean(history.history['mean_absolute_error'])
sum(history.history['mean_squared_error'])/(len(history.history['mean_squared_error']))
nn = history
nncolumn = [ np.mean(nn.history['mean_absolute_error']), np.mean(nn.history['mean_squared_error']), "not calculated", np.mean(nn.history['mean_absolute_percentage_error'])]

eNet = sklearn.linear_model.ElasticNet()
models = [mlr_reg, rf_reg, rf_tuned, xgb_reg, xgb_tuned, lasso_reg, lassotuned, ridge_reg, ridgetuned, eNet, eNettuned, nn]
model_names = ['Linear Regression', 'Random Forest', 'Random Forest (Tuned)', "XGBoost", "XGBoost Tuned", 'Lasso', 'Lasso (Tuned)', 'Ridge', 'Ridge (Tuned)', 'Elastic Net', 'Elastic Net (Tuned)']

# Run the Cross-Validation comparison with the models used in this analysis
comp, maes, mses, r2s, accs = cv_comparison(models, concreteX, concreteY, 5, model_names)
from tabulate import tabulate

print (tabulate(comp, headers= model_names))

# Create DataFrame for all R^2s
r2_comp = pd.DataFrame(r2s, index=comp.columns, columns=['1st Fold', '2nd Fold', '3rd Fold', 
                                                         '4th Fold', '5th Fold'])

# Add a column for the averages
r2_comp['Average'] = np.round(r2_comp.mean(axis=1),4)
r2_comp

# Create DataFrame for all R^2s
df_comp = pd.DataFrame(comp, index=comp.index, columns=model_names)
df_comp["Neural Network"] = nncolumn

# Add a column for the averages
df_comp<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11922268907563026, 'fraction_numerical': 0.025577731092436974, 'mean_word_length': 4.2425660792951545, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 57, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '206722', 'license': 'MIT', 'n_tokens_mistral': 7251, 'n_tokens_neox': 6824, 'n_words': 1585, 'path': 'Lab3_ipnyb.ipynb', 'repo_name': 'NewFigton/newfigton.github.io', 'size': '19040'}