{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import pandas as pd\nbuild_meta = pd.read_csv(\"../input/ashrae-energy-prediction/building_metadata.csv\")\nsample_submission = pd.read_csv(\"../input/ashrae-energy-prediction/sample_submission.csv\")\ntest = pd.read_csv(\"../input/ashrae-energy-prediction/test.csv\")\ntrain = pd.read_csv(\"../input/ashrae-energy-prediction/train.csv\")\nweather_test = pd.read_csv(\"../input/ashrae-energy-prediction/weather_test.csv\")\nweath_train = pd.read_csv(\"../input/ashrae-energy-prediction/weather_train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install rfpimp","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"###########Import required packages###########\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport csv as csv\nimport xgboost as xgb\nfrom xgboost import plot_importance\nfrom matplotlib import pyplot\nfrom sklearn.model_selection import  train_test_split, RandomizedSearchCV\nfrom sklearn.metrics import median_absolute_error, mean_squared_error, mean_absolute_error\nfrom sklearn.metrics import accuracy_score\nfrom scipy.stats import uniform, randint\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.compose import ColumnTransformer\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.preprocessing import LabelEncoder\nfrom rfpimp import *\nfrom sklearn.tree import export_graphviz\nfrom subprocess import call\nfrom IPython.display import Image\nsns.set()\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"###########Data Preprocessing###########\ntrain.timestamp = pd.to_datetime(train.timestamp)\nweath_train.timestamp = pd.to_datetime(weath_train.timestamp)\n\nweath_train['month'] = weath_train['timestamp'].dt.month\nweath_train['day'] = weath_train['timestamp'].dt.day\nweath_train_togrp = weath_train.drop(['timestamp'], axis=1)\nweath_train_daily = weath_train_togrp.groupby(['site_id','month','day']).mean().reset_index()\n\ntrain['day_of_week'] = train['timestamp'].dt.dayofweek\ntrain['month'] = train['timestamp'].dt.month\ntrain['day'] = train['timestamp'].dt.day\ntrain_togrp = train.drop(['timestamp'], axis=1)\ntrain_daily = train_togrp.groupby(['building_id','meter','month','day']).mean().reset_index()\n\nbuild_train_merged = pd.merge(build_meta, train_daily, on='building_id', how='inner')\ndata = pd.merge(build_train_merged, weath_train_daily, on=['site_id','month','day'], how='inner')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"* Here, we identify columns with nulls. We will not use the year_built and floor_count variables because we are hesitant to impute data for variables where majority of the data is missing."},{"metadata":{"trusted":true},"cell_type":"code","source":"# To count the NULL/NaN values and drop columns \nlen(data) - data.count().sort_values(ascending=True)\npercent_missing = data.isnull().sum() * 100 / len(data)\nprint(percent_missing.sort_values(ascending=False))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"########### Split Dataset into Train and Test###########\nX_train, X_test, y_train, y_test = train_test_split(data.drop(columns=['meter_reading']), \n                                                    data[['meter_reading']], \n                                                    test_size=0.25, \n                                                    random_state=42, shuffle=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.shape, y_train.shape, X_test.shape, y_test.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"###########Data Visualization###########\nfeatures = ['site_id','building_id','square_feet','meter','month','air_temperature','dew_temperature',\n            'wind_speed','cloud_coverage','wind_direction','day','precip_depth_1_hr']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Correlation matrix\ndef plotCorrelationMatrix(df, graphWidth):\n#    filename = df.dataframeName\n    df = df.dropna('columns') # drop columns with NaN\n    df = df[[col for col in df if df[col].nunique() > 1]] # keep columns where there are more than 1 unique values\n    if df.shape[1] < 2:\n        print(f'No correlation plots shown: The number of non-NaN or constant columns ({df.shape[1]}) is less than 2')\n        return\n    corr = df.corr()\n    plt.figure(num=None, figsize=(graphWidth, graphWidth), dpi=80, facecolor='w', edgecolor='k')\n    corrMat = plt.matshow(corr, fignum = 1)\n    plt.xticks(range(len(corr.columns)), corr.columns, rotation=90)\n    plt.yticks(range(len(corr.columns)), corr.columns)\n    plt.gca().xaxis.tick_bottom()\n    plt.colorbar(corrMat)\n    plt.title(f'Correlation Matrix for train dataset', fontsize=15)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plotCorrelationMatrix(X_train,6)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(pd.concat([X_train,y_train], axis=1), x_vars=features[:4], y_vars='meter_reading')\nsns.pairplot(pd.concat([X_train,y_train], axis=1), x_vars=features[4:8], y_vars='meter_reading')\nsns.pairplot(pd.concat([X_train,y_train], axis=1), x_vars=features[8:], y_vars='meter_reading')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> We will use label encoding technique for the categorical feature primary_use."},{"metadata":{"trusted":true},"cell_type":"code","source":"le = LabelEncoder()\nX_train.primary_use = le.fit_transform(X_train['primary_use'])\nX_test.primary_use = le.transform(X_test['primary_use'])\nX_train.columns = [col.rstrip('_') for col in X_train.columns] \nX_test.columns = [col.rstrip('_') for col in X_test.columns]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"3. XgBoost\n\nThis histogram of the target variable clearly shows that its long tailed. Squared error would not be a good loss function."},{"metadata":{},"cell_type":"markdown","source":"We want to optimize for the median of the errors since the distribution is long tailed. We will use a Pseudo Huber loss function (shown below), instead of MSE which is very sensitive to outliers.\n\n$L_h(d) = h^2 (\\sqrt{1 + (d/h)^2} - 1)$ where $h = 1$\n\nOur implementation was inspired by this stackoverflow post https://stackoverflow.com/questions/45006341/xgboost-how-to-use-mae-as-objective-function, https://www.kaggle.com/stuarthallows/using-xgboost-with-scikit-learn and a github project by samarthinani"},{"metadata":{"trusted":true},"cell_type":"code","source":"def huber_approx_obj(train, preds):\n    \"\"\"\n    Function returns gradient and hessein of the Pseudo-Huber function.\n    \"\"\"\n    d = preds - train\n    h = 1  ## constant\n    scale = 1 + (d / h) ** 2\n    scale_sqrt = np.sqrt(scale)\n    grad = d / scale_sqrt\n    hess = 1 / scale / scale_sqrt\n    return grad, hess","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## define huber loss - minimizing it means maximizing its negative\ndef huber_loss(preds, train):\n    \"\"\"Function returns the huber loss for h = 1\"\"\"\n    d = preds - train\n    h = 1\n    return -1 * np.sum(np.sqrt(1 + (d/h)**2) - 1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"XgBoost is a very powerful algorithm when it comes to tabular data. Hence we decided to try out the XGBRegressor model.\n\nTo find out the optimal hyperparameter combination , we will do a randomized search over a hyper-dimensional space, fitting XgBoost models to minimize the validation Pseudo-Huber loss.\n\nWe will be fitting a total of 50 XgBoost models by performing a 5-fold cross validation on the train dataset. This will help us improve the generality of the model.\n\nLets wrap up all this functionality into a pipeline and conduct a randomized search."},{"metadata":{"trusted":true},"cell_type":"code","source":"num_features_xgb = ['site_id', 'building_id', 'primary_use', 'square_feet', 'meter', 'month',\n                'day', 'air_temperature', 'dew_temperature', 'wind_speed', 'wind_direction',\n                'cloud_coverage', 'precip_depth_1_hr']\n\nnum_transformer_xgb = Pipeline(steps=[('imputer', SimpleImputer(strategy='median'))])\n\npreprocessor_xgb = ColumnTransformer(transformers=[('num', num_transformer_xgb, num_features_xgb)])\n\nparams_xgb = {\n    \"rcv_xgb__colsample_bytree\": uniform(0.7, 0.1),\n    \"rcv_xgb__gamma\": uniform(0, 0.2),\n    \"rcv_xgb__learning_rate\": uniform(0.03, 0.12), \n    \"rcv_xgb__subsample\": uniform(0.8, 0.15),\n    \"rcv_xgb__booster\": ['gbtree','dart']\n}\n\npipeline_xgb = Pipeline(steps=[('preprocessor', preprocessor_xgb),\n                           ('rcv_xgb', xgb.XGBRegressor(objective=huber_approx_obj, \n                                                        feval= huber_loss, max_depth=5,n_estimators=30))])\nsearch_xgb = RandomizedSearchCV(pipeline_xgb, param_distributions=params_xgb, n_iter=10, \n                            scoring='neg_median_absolute_error', random_state=42, cv=5, \n                            verbose=1, n_jobs=4, return_train_score=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nsearch_xgb.fit(X_train[num_features_xgb], y_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Evaluation Metrics**"},{"metadata":{"trusted":true},"cell_type":"code","source":"def report_best_scores(results, n_top=3):\n    \"\"\"Function gives hyperparameters for the top n models\"\"\"\n    for i in range(1, n_top + 1):\n        candidates = np.flatnonzero(results['rank_test_score'] == i)\n        for candidate in candidates:\n            print(\"Model with rank: {0}\".format(i))\n            print(\"Mean validation score: {0:.3f} (std: {1:.3f})\".format(\n                  results['mean_test_score'][candidate],\n                  results['std_test_score'][candidate]))\n            print(\"Parameters: {0}\".format(results['params'][candidate]))\n            print(\"\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"report_best_scores(search_xgb.cv_results_, 3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"xgb_predictions = search_xgb.predict(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train.hist(bins = 100, range = [1,2000]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Evaluation Metric**\n* A histogram of the target variable is long tailed. We will use the Median Absolute Error as our North Star evaluation metric."},{"metadata":{"trusted":true},"cell_type":"code","source":"y_train_pred_xgb = search_xgb.predict(X_train[num_features_xgb]) \ny_test_pred_xgb = search_xgb.predict(X_test[num_features_xgb])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_medae_xgb = median_absolute_error(y_train, y_train_pred_xgb)\ntest_medae_xgb = median_absolute_error(y_test, y_test_pred_xgb)\nprint(f'MEDAE XGBoost: Train = {train_medae_xgb:.2f} , Test = {test_medae_xgb:.2f}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The XGBoost model, being the most powerful for tabular data, gave us the best test MedAE value of 48.\n\nWe also looked at the RMSE and MAE for our reference."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_rmse_xgb = np.sqrt(mean_squared_error(y_train, y_train_pred_xgb))\ntest_rmse_xgb = np.sqrt(mean_squared_error(y_test, y_test_pred_xgb))\nprint(f'RMSE XGBoost: Train = {train_rmse_xgb:.2f} , Test = {test_rmse_xgb:.2f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_mae_xgb = mean_absolute_error(y_train, y_train_pred_xgb)\ntest_mae_xgb = mean_absolute_error(y_test, y_test_pred_xgb)\nprint(f'MAE XGBoost: Train = {train_mae_xgb:.2f} , Test = {test_mae_xgb:.2f}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_medae = huber_loss(y_train.values.ravel(), y_train_pred_xgb)\ntest_medae = huber_loss(y_test.values.ravel(), y_test_pred_xgb)\nprint(f'Huber Loss XGBoost: Train = {train_medae:.2f} , Test = {test_medae:.2f}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Feature Importances\n![](http://)**\nThe square feet area of the buidling and its meter type turn out to be two the most important variables"},{"metadata":{"trusted":true},"cell_type":"code","source":"I_xgb = importances(search_xgb.best_estimator_, X_test[num_features_xgb], y_test)\nprint(I_xgb)\nplot_importances(I_xgb,title= 'Feature Importance',imp_range=(0, 0.05))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test.reset_index(inplace=True)\nX_test['meter_reading'] = xgb_predictions[:,]\nsample_submission = X_test[['index','meter_reading']]\nsample_submission = sample_submission.rename(columns={\"index\": \"row_id\"})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"########Predicted Results#########\nprint(sample_submission)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission.to_csv('sample_submission.csv')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}