CO Prediction Code.py

# -*- coding: utf-8 -*-
"""ML Assignment Code - Geoffrey Hinton .ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1IZMc3fvrvhWYKYZf5hj4DhE3EfSA5B-z

**Boruta installation -** Boruta needs to be installed in order to do the feature selection using BorutaPy
"""

pip install boruta

"""**Import all necessary libraries**"""

# Libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime
from sklearn import preprocessing
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import GridSearchCV
from sklearn import metrics
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import RandomForestClassifier
from sklearn.preprocessing import StandardScaler
from sklearn import decomposition, datasets
from sklearn import linear_model
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import LabelEncoder
from boruta import BorutaPy

#Read the dataset
df = pd.read_csv("city_day.csv")
df

"""**Exploratory Data Analysis**"""

print ('Rows with no. of missing values:')
print()
print(df.isna().sum())

#Information about the dataframe
df.info()

df.describe()

#Outliers
sns.boxplot(data=df[['PM2.5','PM10','NO','NO2','NOx','NH3','CO','SO2','O3','Benzene','Toluene','Xylene','AQI' ]], orient='h')
plt.show()

sns.heatmap(df.corr(), annot = True, cmap = 'magma')

"""**Data Cleaning**"""

#Dropping attribute that do not contribute to the findings
df.drop(['AQI_Bucket'], axis=1, inplace = True)
df

# #Drop the rows with null values as the number of missing values not large enough to affect the size of the dataset
df2 = df.dropna()
df2.head()

#convert the values to 32-bit integer data type
df2['PM2.5'] = df2['PM2.5'].apply(np.int32)
df2['PM10'] = df2['PM10'].apply(np.int32)
df2['NO'] = df2['NO'].apply(np.int32)
df2['NO2'] = df2['NO2'].apply(np.int32)
df2['NOx'] = df2['NOx'].apply(np.int32)
df2['NH3'] = df2['NH3'].apply(np.int32)
df2['CO'] = df2['CO'].apply(np.int32)
df2['SO2'] = df2['SO2'].apply(np.int32)
df2['O3'] = df2['O3'].apply(np.int32)
df2['Benzene'] = df2['Benzene'].apply(np.int32)
df2['Toluene'] = df2['Toluene'].apply(np.int32)
df2['Xylene'] = df2['Xylene'].apply(np.int32)
df2['AQI'] = df2['AQI'].apply(np.int32)

#Label Encoding done to attribute 'City'
le = LabelEncoder()
df2['City'] = le.fit_transform(df2['City'])
df2

y = df2['CO']
X = df2.drop(['CO','Date'] , 1)
colnames = X.columns

def ranking(ranks, names, order=1):
    minmax = MinMaxScaler()
    ranks = minmax.fit_transform(order*np.array([ranks]).T).T[0]
    ranks = map(lambda x: round(x,2), ranks)
    return dict(zip(names, ranks))

rf = RandomForestClassifier(n_jobs=-1,class_weight="balanced",max_depth = 5)

feat_selector=BorutaPy(rf,n_estimators="auto",random_state = 1)

feat_selector.fit(X.values,y.values.ravel())
# sort the boruta score descending

boruta_score = ranking(list(map(float, feat_selector.ranking_)), colnames, order=-1)
boruta_score = pd.DataFrame(list(boruta_score.items()), columns=['Features', 'Score'])
boruta_score = boruta_score.sort_values("Score",ascending= False)

print('---------Top 10----------')
# your codes here...
display(boruta_score.head(10))

print('---------Bottom 10----------')
# your codes here...
boruta_score.tail(10)

predictors = ['PM2.5','PM10','NO','NOx','AQI','NH3','Benzene','SO2','O3','NO2','Toluene','Xylene']
X = df2[predictors]
y = df2['CO']
colnames = X.columns

# Splitting the dataset into training and testing set (70/30)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.3, random_state = 0)

"""#**Decision Tree Regression**"""

#imports packages used for making decision tree and tuning its hyperparameters
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import GridSearchCV, train_test_split
from sklearn.metrics import r2_score

UntunedDecisionTree=DecisionTreeRegressor()
UntunedDecisionTree.fit(X_train,y_train)
dtr_prediction = UntunedDecisionTree.predict(X_test)

print('MAE Test:', metrics.mean_absolute_error(y_test,dtr_prediction))
print('RMSE Test:', np.sqrt(metrics.mean_squared_error(y_test, dtr_prediction)))
print('R^2:',r2_score(y_test, dtr_prediction))

r2_dt= r2_score(y_test, dtr_prediction)
N = len(y_test)
k = 3
r2_dts_3 = (1-r2_dt)*(N-1)/(N-k-1)
print('Adjusted R-Squared =',r2_dts_3)

r2_dt= r2_score(y_test, dtr_prediction)
N = len(y_test)
k = 6
r2_dts_6 = (1-r2_dt)*(N-1)/(N-k-1)
print('Adjusted R-Squared =',r2_dts_6)

r2_dt= r2_score(y_test, dtr_prediction)
N = len(y_test)
k = 9
r2_dts_9 = (1-r2_dt)*(N-1)/(N-k-1)
print('Adjusted R-Squared =',r2_dts_9)

"""### Tuning the Model by optimizing Hyperparameters utilizing GridSearchCV"""

#Grid of parameters of DecisionTreeRegression() and their respective values
#Parameters with input values of type Int or Float have values chosen randomly to be put in the grid. Value is then updated/replaced based on optimal value of last iteration
#Parameters with a small fixed set of values such as 'max_features' and 'criterion' have all values included in the grid
#https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html - List of parameters for DecisionTreeRegressor() function

param_grid = {'max_depth': [23,24,25,26],                                  #Max Depth of 26 as data originally had 26 features before being cut
              'min_samples_split': [14,15,16,17],                          #min samples for decision node to avoid overfitting
              'min_samples_leaf': [5,6,8,10],                              #min amount of samples in leaf node to avoid overfitting
              'max_features':[2,5,8],                                      #max number of unique features for leaf node
              'criterion':["squared_error","friedman_mse", "entropy"],
              'min_impurity_decrease': [0.2, 0.4, 0.6, 0.8]}               #min impurity measure to maximise information gain

#Tuning model so that the decision tree doesnt overfit or underfit
#Uses the param_grid defined and iterates through every single combination to find best possible
#https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html

DTR_TunedModel = GridSearchCV(DecisionTreeRegressor(),
                   param_grid,
                   cv=3, n_jobs=-1, verbose=2, pre_dispatch='2*n_jobs')
DTR_TunedModel.fit(X_train,y_train)

DTR_TunedModel.best_estimator_, DTR_TunedModel.best_params_, DTR_TunedModel.best_score_

dtr_tuned_prediction = DTR_TunedModel.predict(X_test)
plt.scatter(y_test, dtr_tuned_prediction)

#plot of untuned model prediction for comparison
plt.scatter(y_test, dtr_prediction)

print('MAE Test:', metrics.mean_absolute_error(y_test,dtr_tuned_prediction))
print('RMSE Test:', np.sqrt(metrics.mean_squared_error(y_test,dtr_tuned_prediction)))
print('R^2:',r2_score(y_test,dtr_tuned_prediction))

"""# **Support Vector Regression**"""

from sklearn.metrics import mean_squared_error
import math

from sklearn.svm import SVR

svr = SVR()
svr.fit(X_train, y_train)

y_pred = svr.predict(X_test)

from sklearn.metrics import r2_score
print('MAE Test:', metrics.mean_absolute_error(y_test,y_pred))
print('RMSE Test:', np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
print('R^2:',r2_score(y_test, y_pred))

from sklearn.metrics import r2_score
r2_svr = r2_score(y_test, y_pred)
N = len(y_test)
k = 3
r2_svr_3 = (1-r2_svr)*(N-1)/(N-k-1)
print('Adjusted R-Squared =', r2_svr_3)

from sklearn.metrics import r2_score
r2_svr = r2_score(y_test, y_pred)
N = len(y_test)
k = 6
r2_svr_6= (1-r2_svr)*(N-1)/(N-k-1)
print('Adjusted R-Squared =', r2_svr_6)

r2_svr = r2_score(y_test, y_pred)
N = len(y_test)
k = 9
r2_svr_9 = (1-r2_svr)*(N-1)/(N-k-1)
print('Adjusted R-Squared =', r2_svr_9)

"""### Tuning the Model by optimizing Hyperparameters utilizing GridSearchCV"""

from sklearn.linear_model import Ridge
from sklearn.model_selection import RepeatedKFold
from numpy import arange

# define model
model = Ridge()
# define model evaluation method
cv = RepeatedKFold(n_splits=10, n_repeats=3, random_state=1)

# define grid
grid = dict()
grid['alpha'] = arange(0, 1, 0.01)
# define search
search = GridSearchCV(model, grid, scoring='neg_mean_absolute_error', cv=cv, n_jobs=-1)
# perform the search
search.fit(X_train,y_train)

search.best_params_

predictions_MLR = search.predict(X_test)

print('MAE Test:', metrics.mean_absolute_error(y_test,predictions_MLR))
print('RMSE Test:', np.sqrt(metrics.mean_squared_error(y_test,predictions_MLR)))
print('R^2:',r2_score(y_test,predictions_MLR))

"""# **Deep Learning Model - Convolutional Neural Networks (CNN)**"""

pip install keras-tuner -q

#import required TensorFlow libraries
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from sklearn.preprocessing import StandardScaler
from kerastuner.tuners import RandomSearch
from kerastuner.engine.hyperparameters import HyperParameter

"""**Building the model and tuning the  Hyperparameters using KerasTuner**"""

#Standardize the input features
scaler = StandardScaler()

X = scaler.fit_transform(df2[predictors])
y = df2['CO']

#Define the model
def build_model(hp):
  model = Sequential()

  #Tuning the no. of units in the 2 dense layers
  hp_units = hp.Int('units', min_value = 32, max_value = 512, step = 32)
  model.add(Dense(units = hp_units, activation = 'relu'))

  hp_units = hp.Int('units', min_value = 32, max_value = 256, step = 32)
  model.add(Dense(units = hp_units, activation = 'relu'))

  #Output layer with linear activation function for regression
  model.add(Dense(1, activation = 'linear'))

  #Tuning the learning rate for the optimizer
  hp_learning_rate = hp.Choice('learning_rate', values = [1e-2, 1e-3, 1e-4])

  model.compile(loss = 'mean_absolute_error', optimizer = keras.optimizers.Adam(learning_rate = hp_learning_rate))

  return model

#Initializing the tuner using RandomSearch
tuner = RandomSearch(
    build_model,
    objective = 'val_loss',
    max_trials = 10,
    directory = 'tuner_dir',
    project_name = 'co_tuner'
)

#Hyperparameter search
tuner.search(X_train, y_train, epochs = 10, batch_size = 32, validation_split = 0.3)

#Obtain the best hyperparameters
best_hp = tuner.get_best_hyperparameters(num_trials = 1)[0]

#Final model with best hyperparameters
final_model = build_model(best_hp)

final_model.fit(X_train, y_train, epochs = 10, batch_size = 32, validation_split = 0.3)

#Evaluate final model performance using MAE, RMSE & R^2 values

#Prediction of target values using the final model
y_pred = final_model.predict(X_test)

mae = final_model.evaluate(X_test, y_test)
print(f'Test MAE: {mae}')

#RMSE Calculation
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
print(f'Test RMSE: {rmse}')

#R^2 Value Calculation
r2 = r2_score(y_test, y_pred)
print(f'Test R^2: {r2}')