<jupyter_start><jupyter_text>============== Muti-layer Perceptron (publically available dataset)  ===============st8_Zhang_Dong_codeCreated on 1st July, Thu, 2021- dataset: energydata_complete.csv (See the report for dataset details)- To predict the energy consumption- regression problem please read the readme.txt file first.<jupyter_code># import packages
from sklearn import datasets
from sklearn import metrics
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import cross_val_score
from sklearn.exceptions import ConvergenceWarning
from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
# always print floating point numbers using fixed point notation
np.set_printoptions(suppress=True)

# set plotting style
plt.style.use('ggplot')

# import warnings filter
from warnings import simplefilter

# ignore all future warnings and Convergence Warning
simplefilter(action='ignore', category=FutureWarning)
simplefilter(action='ignore', category=ConvergenceWarning)
# load data
# set the first row as head and the first column as index
df = pd.read_csv('../energydata_complete.csv', header=0, index_col=0)
df.head()
# convert a pandas Dataframe to numpy array
df = np.array(df)
# get the length of dataset
nums = df.shape[1]
# display the demo of dataset
print('The first three rows of the dataset:\n', df[:3,:])
print('\nThe total columns of dataset is: {}'.format(nums))
# Modify df in-place by shuffling its contents.
np.random.shuffle(df)

# split the dataset: train_set & test_set
ratio = 0.8
offset = int(df.shape[0] * ratio)

train_data = df[:offset]

# normalize data
std = train_data.std(axis=0)
mean = train_data.mean(axis=0)

for i in range(nums):
    # x <- (x-u)/s
    df[:, i] = (df[:, i] - mean[i]) / (std[i])

train_data = df[:offset]
test_data = df[offset:]

# make the first column "crim" as the target
X_train = train_data[:, 1:]  # training set: data
y_train = train_data[:, 0]   # training set: target
X_test = test_data[:, 1:]    # testing set: data
y_test = test_data[:, 0]     # testing set: target
print('The shape of X_train: ', X_train.shape)
print(X_train[0])
print('The shape of y_train: ', y_train.shape)
print(y_train[0])
print('The shape of X_test ', X_test.shape)
print(X_test[0])
print('The shape of  y_test: ', y_test.shape)
print(y_test[0])<jupyter_output>The shape of  y_test:  (3947,)
-0.6596953202694181<jupyter_text>===================== Multilayer perceptron ======================<jupyter_code># MLP model
# alpha_l1: the alpha for L1 regularization
# alpha_l2: the alpha for L2 regularization
# to set alpha_l1=0 and alpha_l2=0, so there is no regularization used in the model
model_regression = MLPRegressor(alpha_l1 = 0.0, alpha_l2=0.0)
# Train the model
model_regression.fit(X_train, y_train)

# train set
expected_x = y_train
predicted_x = model_regression.predict(X_train)
print('MLP: The R2 score for training set is: ', metrics.r2_score(expected_x, predicted_x))
print('MLP: MSE on training set:', metrics.mean_squared_error(expected_x, predicted_x))
print('')

# test set
expected_y = y_test
predicted_y = model_regression.predict(X_test)
print('MLP: The R2 score for testing is: ', metrics.r2_score(expected_y, predicted_y))
print('MLP: MSE on testing set:', metrics.mean_squared_error(expected_y, predicted_y))<jupyter_output>MLP: The R2 score for training set is:  0.5027987998932935
MLP: MSE on training set: 0.4972012001067066

MLP: The R2 score for testing is:  0.35995641751569174
MLP: MSE on testing set: 0.6197757805003244<jupyter_text>================== Multilayer perceptron with Lasso ===================<jupyter_code># ================ MLP with L1 regularization ==================
# train alpha_l1
# create an empty list to store training result
l1_param_train = [[], [], []]   # alpha_l1, score_l1_R2, score_l1_RME
l1_param_test = [[], [], []]   # alpha_l1, score_l1_R2, score_l1_RME

expected_l1_x = y_train
expected_l1_y = y_test

# Here, the idea of grid search is applied to train Lasso model 
# to find the optimal alpha value in the interval of (0,2]
# step length is 0.1
for alpha_train_l1 in range(1,21):
    # using Grid Search to search the best alpha in [0,2]
    # MLP model with L1 penalty
    # in lasso model, set alpha_l2 = 0.0
    mlp_l1 = MLPRegressor(alpha_l1 = alpha_train_l1*0.1, alpha_l2=0.0)
    mlp_l1.fit(X_train, y_train)
    
    # training set
    predicted_l1_x = mlp_l1.predict(X_train)
    s_l1_train_r2 = metrics.r2_score(expected_l1_x, predicted_l1_x)
    s_l1_train_mse = metrics.mean_squared_error(expected_l1_x, predicted_l1_x)
    # append training results to the l1_param_train
    l1_param_train[0].append(alpha_train_l1*0.1)
    l1_param_train[1].append(s_l1_train_r2)
    l1_param_train[2].append(s_l1_train_mse)
    
    # testing set
    predicted_l1_y = mlp_l1.predict(X_test)
    s_l1_test_r2 = metrics.r2_score(expected_l1_y, predicted_l1_y)
    s_l1_test_mse = metrics.mean_squared_error(expected_l1_y, predicted_l1_y)
    # append testing results to the l1_param_test
    l1_param_test[0].append(alpha_train_l1*0.1)
    l1_param_test[1].append(s_l1_test_r2)
    l1_param_test[2].append(s_l1_test_mse)
    
%matplotlib inline
# plot alpha vs. R^2
plt.title('MLP with L1')
plt.plot(l1_param_train[0], l1_param_train[1], label='R^2_train')
plt.plot(l1_param_train[0], l1_param_test[1],  label='R^2_test')
plt.xlabel("alpha")
plt.legend()
plt.show()
%matplotlib inline
# plot alpha vs. RMSE
plt.title('MLP with L1')
plt.plot(l1_param_train[0], l1_param_train[2], label='MSE_train')
plt.plot(l1_param_train[0], l1_param_test[2],  label='MSE_test')
plt.xlabel("alpha")
plt.legend()
plt.show()
# print the result of MLP_l1 on testing data in (0,2]
l1_param_test = np.array(l1_param_test)
# get the index of minimun mse, R^2 = (1-SSR)/SST, so R^2 is maximum when MSE is minimum
a_l1_test = np.argmin(l1_param_test, axis=1)
# print(a_l1_test)
print('The best result for MLP with L1 regularization: \nalpha = {} \nR^2 = {} \nMSE = {}'
      .format(l1_param_test[0][a_l1_test[2]], l1_param_test[1][a_l1_test[2]], l1_param_test[2][a_l1_test[2]]))<jupyter_output>The best result for MLP with L1 regularization: 
alpha = 0.1 
R^2 = 0.36250160191781766 
MSE = 0.6173111926308202<jupyter_text>====================== MLP with Ridge ======================<jupyter_code># train alpha_l2
# create an empty list to store training result
l2_param_train = [[], [], []]   # alpha_l2, score_l2_R2, score_l2_RME
l2_param_test = [[], [], []]   # alpha_l2, score_l2_R2, score_l2_RME

expected_l2_x = y_train
expected_l2_y = y_test

# grid search
# to find the optimal alpha value in the interval of (0,2]
# step length is 0.1
for alpha_train_l2 in range(1,21):
    # using Grid Search to search the best alpha in [0,2]
    # MLP model with L1 penalty
    # in ridge model, set alpha_l1 = 0.0
    mlp_l2 = MLPRegressor(alpha_l1 = 0.0, alpha_l2=alpha_train_l2*0.1)
    mlp_l2.fit(X_train, y_train)
    
    # training set
    predicted_l2_x = mlp_l2.predict(X_train)
    s_l2_train_r2 = metrics.r2_score(expected_l2_x, predicted_l2_x)
    s_l2_train_mse = metrics.mean_squared_error(expected_l2_x, predicted_l2_x)
    # append training results to the l2_param_train
    l2_param_train[0].append(alpha_train_l2*0.1)
    l2_param_train[1].append(s_l2_train_r2)
    l2_param_train[2].append(s_l2_train_mse)
    
    # testing set
    predicted_l2_y = mlp_l2.predict(X_test)
    s_l2_test_r2 = metrics.r2_score(expected_l2_y, predicted_l2_y)
    s_l2_test_mse = metrics.mean_squared_error(expected_l2_y, predicted_l2_y)
    # append testing results to the l2_param_train
    l2_param_test[0].append(alpha_train_l2*0.1)
    l2_param_test[1].append(s_l2_test_r2)
    l2_param_test[2].append(s_l2_test_mse)
%matplotlib inline
# plot alpha vs. R^2
plt.title('MLP with L2')
plt.plot(l2_param_train[0], l2_param_train[1], label='R^2_train')
plt.plot(l2_param_train[0], l2_param_test[1],  label='R^2_test')
plt.xlabel("alpha")
plt.legend()
plt.show()
%matplotlib inline
# plot alpha vs. RMSE
plt.title('MLP with L2')
plt.plot(l2_param_train[0], l2_param_train[2], label='MSE_train')
plt.plot(l2_param_train[0], l2_param_test[2],  label='MSE_test')
plt.xlabel("alpha")
plt.legend()
plt.show()
# print the result of MLP_l2 on testing data in [0,2]
l2_param_test = np.array(l2_param_test)
# get the index of minimun mse
a_l2_test = np.argmin(l2_param_test, axis=1)
# print(a_l2_test)
print('The best result for MLP with L2 regularization: \nalpha = {} \nR^2 = {} \nMSE = {}'
      .format(l2_param_test[0][a_l2_test[2]], l2_param_test[1][a_l2_test[2]], l2_param_test[2][a_l2_test[2]]))<jupyter_output>The best result for MLP with L2 regularization: 
alpha = 0.1 
R^2 = 0.37656001785207827 
MSE = 0.60369795449095<jupyter_text>============== MLP with elastic net ================<jupyter_code># train alpha_l1_l2
# create an empty list to store training result
l1_l2_param_train = [[], [], [], []]   # alpha_l1_l2, score_l1_l2_R2, score_l1_l2_RME
l1_l2_param_test = [[], [], [], []]   # alpha_l1_l2, score_l1_l2_R2, score_l1_l2_RME

expected_l1_l2_x = y_train
expected_l1_l2_y = y_test

# grid search
# to find the optimal alpha value in the interval of (0,2]
# step length is 0.1

for alpha_l1 in range(1,21):
    for alpha_l2 in range(1,21):
        # using Grid Search to search the best alpha in [0,2]
        # MLP model with L1 & L2 penalty
        mlp_l1_l2 = MLPRegressor(alpha_l1=alpha_l1*0.1, alpha_l2=alpha_l2*0.1)
        mlp_l1_l2.fit(X_train, y_train)
        
        # training set
        predicted_l1_l2_x = mlp_l1_l2.predict(X_train)
        s_l1_l2_train_r2 = metrics.r2_score(expected_l1_l2_x, predicted_l1_l2_x)
        s_l1_l2_train_mse = metrics.mean_squared_error(expected_l1_l2_x, predicted_l1_l2_x)
        # append training results to the l1_l2_param_train
        l1_l2_param_train[0].append(alpha_l1*0.1)
        l1_l2_param_train[1].append(alpha_l2*0.1)
        l1_l2_param_train[2].append(s_l1_l2_train_r2)
        l1_l2_param_train[3].append(s_l1_l2_train_mse)

        # testing set
        predicted_l1_l2_y = mlp_l1_l2.predict(X_test)
        s_l1_l2_test_r2 = metrics.r2_score(expected_l1_l2_y, predicted_l1_l2_y)
        s_l1_l2_test_mse = metrics.mean_squared_error(expected_l1_l2_y, predicted_l1_l2_y)
        # append testing results to the l1_l2_param_test
        l1_l2_param_test[0].append(alpha_l1*0.1)
        l1_l2_param_test[1].append(alpha_l2*0.1)
        l1_l2_param_test[2].append(s_l1_l2_test_r2)
        l1_l2_param_test[3].append(s_l1_l2_test_mse)
# you can save the result to a local .csv file by the following code
# np.savetxt("MLP_L1&L2_energy.csv", l1_l2_param_test.T, delimiter =",",fmt ='% s')
# print the result of MLP_l1_l2 on testing data in [0,2]
l1_l2_param_test = np.array(l1_l2_param_test)
a_l1_l2_test = np.argmin(l1_l2_param_test, axis=1)
# print(a_l1_l2_test)
print('The best result for MLP with L1&L2 regularization: \nalpha_l1 = {} \nalpha_l2 = {} \nR^2 = {} \nMSE = {}'
      .format(l1_l2_param_test[0][a_l1_l2_test[3]], l1_l2_param_test[1][a_l1_l2_test[3]], 
              l1_l2_param_test[2][a_l1_l2_test[3]], l1_l2_param_test[3][a_l1_l2_test[3]]))
# plot the result of MLP with L1&L2 regularization
%matplotlib inline
fig = plt.figure()
ax = Axes3D(fig)
#  to make "-" show normally
plt.rcParams['axes.unicode_minus']=False

# plotting
surf = ax.scatter(l1_l2_param_test[0],l1_l2_param_test[1],l1_l2_param_test[2],cmap="Blues")
ax.view_init(elev=15, azim=15)

# setting tittle and labels
plt.title('MLP with L1 & L2: testing set R^2')
ax.set_zlabel('score', fontdict={'size': 10, 'color': 'red'})
ax.set_ylabel('alpha_l2', fontdict={'size': 10, 'color': 'red'})
ax.set_xlabel('alpha_l1', fontdict={'size': 10, 'color': 'red'})

# plot the result of MLP with L1&L2 regularization
%matplotlib inline
fig = plt.figure()
ax = Axes3D(fig)

plt.rcParams['axes.unicode_minus']=False #  to make "-" show normally

# plotting
surf = ax.scatter(l1_l2_param_test[0],l1_l2_param_test[1],l1_l2_param_test[3],cmap="Blues")
ax.view_init(elev=25, azim=15)

# setting tittle and labels
plt.title('MLP with L1 & L2: testing set SME')
ax.set_zlabel('score', fontdict={'size': 10, 'color': 'red'})
ax.set_ylabel('alpha_l2', fontdict={'size': 10, 'color': 'red'})
ax.set_xlabel('alpha_l1', fontdict={'size': 10, 'color': 'red'})<jupyter_output><empty_output><jupyter_text>============= Grid search for other parameters ===================<jupyter_code>param_train = [[],[],[],[]]   # alpha_l1, alpha_l2, score_l1_l2_R2, score_l1_l2_RME
param_test = [[],[],[],[]]   # alpha_l1, alpha_l2, score_l1_l2_R2, score_l1_l2_RME

expected_x = y_train
expected_y = y_test

# here, we test two different hidden layer sizes and three diffrent solvers
# hidden_layer_sizes: (100,), (100, 30)
# solver: 'adam', 'sgd', 'lbfgs'
# the value of alpha_l1 and alpha_l2 are from the well-trained model before

for hidden_layer_sizes in [(100,), (100, 30)]:
    for solver in ['adam', 'sgd', 'lbfgs']:
        MLP = MLPRegressor(alpha_l1=1.4, alpha_l2=0.3, 
                                 hidden_layer_sizes=hidden_layer_sizes, solver=solver)
        MLP.fit(X_train, y_train)
        
        # training set
        predicted_x = MLP.predict(X_train)
        s_train_r2 = metrics.r2_score(expected_x, predicted_x)
        s_train_mse = metrics.mean_squared_error(expected_x, predicted_x)
        # append training result to the param_train list
        param_train[0].append(hidden_layer_sizes)
        param_train[1].append(solver)
        param_train[2].append(s_train_r2)
        param_train[3].append(s_train_mse)

        # testing set
        predicted_y = MLP.predict(X_test)
        s_test_r2 = metrics.r2_score(expected_y, predicted_y)
        s_test_mse = metrics.mean_squared_error(expected_y, predicted_y)
        # append testing result to the param_test list
        param_test[0].append(hidden_layer_sizes)
        param_test[1].append(solver)
        param_test[2].append(s_test_r2)
        param_test[3].append(s_test_mse)
# print the result of MLP
param_test = np.array(param_test)
# get the index of minimum mse
a_test = np.argmin(param_test, axis=1)
# print(a_test)
print('The best result for MLP:')
print('alpha_l1 = {} \nalpha_l2 = {}'.format(1.4, 0.3))
print('hidden_layer_sizes = {} \nsolver = {} \nR^2 = {} \nMSE = {}'
      .format(param_test[0][a_test[3]], param_test[1][a_test[3]], 
              param_test[2][a_test[3]], param_test[3][a_test[3]]))
# utilizing the trained parameters
MLP_op = MLPRegressor(alpha_l1=1.4, alpha_l2=0.3, hidden_layer_sizes=(100,30), solver='adam')
MLP_op.fit(X_train, y_train)
expected_op_y = y_test
predicted_op_y = MLP_op.predict(X_test)

# plot the fitting diagram of y and yhat on testing dataset
plt.figure(figsize=(8, 8))
plt.xlabel("y")
plt.ylabel("yhat")
sns.regplot(expected_op_y, predicted_op_y, fit_reg=True, scatter_kws={"s": 30})<jupyter_output><empty_output><jupyter_text>k-fold validation- Below show the k-fold (k=5 & k=10) cross validation for the MLP_op model.- In the previous experiment, the whole data set was divided into training set and test set, so it was no longer necessary to use k-fold for validation. - The results below are only for the purpose of showing the results of model training.<jupyter_code># K-fold cross validation
print('k-fold cross validation:\n')
scores_5 = []
scores_5 = cross_val_score(MLP_op, X_train, y_train, cv=5, scoring='neg_mean_squared_error', n_jobs=1)
print('MSE Loss for k = 5:\n %s\n'% scores_5)

# neg_mean_squared_error return the negative mse, but mse should be a positive number
# fix the sign of MSE scores
mse_scores_5 = -scores_5

print ('Fixed MSE Loss for k = 5:\n %s\n'% mse_scores_5)
print('In conclusion: MSE Loss for k = 5: %3f +/- %3f'% (np.mean(mse_scores_5), np.std(mse_scores_5)))
# K-fold cross validation
print('k-fold cross validation:\n')
scores_10 = []
scores_10 = cross_val_score(MLP_op, X_train, y_train, cv=10, scoring='neg_mean_squared_error', n_jobs=1)
print('MSE Loss for k = 10:\n %s'% scores_10)

# neg_mean_squared_error return the negative mse, but mse should be a positive number
# fix the sign of MSE scores
mse_scores_10 = -scores_10

print ('\nFixed MSE Loss for k = 10:\n %s'% mse_scores_10)
print('\nIn conclusion: MSE Loss for k = 10: %3f +/- %3f'% (np.mean(mse_scores_10), np.std(mse_scores_10)))<jupyter_output>k-fold cross validation:

MSE Loss for k = 10:
 [-0.70599384 -0.67262064 -0.73805399 -0.57447978 -0.58697355 -0.56849987
 -0.72707438 -0.63658375 -0.62574739 -0.66799767]

Fixed MSE Loss for k = 10:
 [0.70599384 0.67262064 0.73805399 0.57447978 0.58697355 0.56849987
 0.72707438 0.63658375 0.62574739 0.66799767]

In conclusion: MSE Loss for k = 10: 0.650402 +/- 0.058949
{'fraction_non_alphanumeric': 0.13195273334924804, 'fraction_numerical': 0.061709238481737885, 'mean_word_length': 4.997494631352899, 'pattern_counts': {'":': 1, '<': 25, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '54852', 'license': 'Apache-2.0', 'n_tokens_mistral': 7114, 'n_tokens_neox': 6687, 'n_words': 1462, 'path': 'data_analytics/MLPRegression_energy.ipynb', 'repo_name': 'CATHY0706/MLP', 'size': '16756'}