{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook is about Kmeans clustering of hourly electricity consumption in Jan 2016. Data is available on Kaggle under Ashrae Energy Prediction Competition. I am interested in finding patterns or perhaps different types of consumers based on hourly consumption.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# import packages.\n\n%matplotlib inline\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom sklearn.cluster import KMeans\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.metrics import silhouette_score\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"## data path\ntrain_path = \"../input/ashrae-energy-prediction/train.csv\"\ndf = pd.read_csv(train_path)\ndf.head","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## usage is in BKTU: convert meter_readings into kwh by multiplying by 0.2931 from BKTU\ndf.loc[:,'meter_reading'] *= 0.2931 ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Data description\n- building_id: type of building grouping different type of users\n- meter: type of meters: 0 for electricity, 1 for hot water\n- timestamp: time at which the reading was registered.\n- meter_reading: usage per hour.\n\nFor the purpose of this analysis, I will only look at one month electricity usage","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## sample data: subset meter == 0 (electricity)\n#train = train.filter(train['meter'] == 0)\ndf = df[(df['meter'] == 0) & (df['timestamp'] <= \"2016-01-31 00:00:00\" )]\n#df.head","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## feature engineering\n ## extract day and month from datestamp column\ndf = df.assign(**{'hour': pd.to_datetime(df['timestamp']).dt.hour,\n                  'day': pd.to_datetime(df['timestamp']).dt.dayofyear,\n                  'day_name': pd.to_datetime(df['timestamp']).dt.day_name()\n                 \n                 }).drop('meter', 1)\n\n#'day': pd.to_datetime(df['timestamp']).dt.day,","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df[(df['meter_reading'] >=8626)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## drop duplicates if there is any\ndf = df.drop_duplicates()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.boxplot(x=df['meter_reading'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The boxplot displays values that can be considered as outliers. However, they might be genuine usage from same building. Let's check that hypothesis by checking if there are one off values (perhaps one hour of the day) or consisent values across the whole day","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# let assume one building can consume 500 kwh.\nhigh_val = df[(df['meter_reading'] >= 500)]\nhigh_val.head","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"high_val['building_id'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# plot high_val: big electricity users.\nhigh_chart = sns.lineplot(x=\"hour\",\n                         y=\"meter_reading\",\n                         data=high_val\n                         ).set_title('hourly distribution for usage above 500')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"high_values data has 35 unique building_ids. Looking at the plot, those values are not outliers, they are spread thoughout the day as normal usage; those buildings might be companies with constant processes such as factories?","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"let's split data into weekends and weekdays: usage patterns might be different.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# split df into weekdays and weekends\n## weekdays: df1\ndf1 = df[(df['day_name'] != 'Saturday') & (df['day_name'] != 'Sunday')]\n\n# weekends: df2\ndf2 = df[(df['day_name'] == 'Saturday') | (df['day_name'] == 'Sunday')]\n# df2.head()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df1.head","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df2.head","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## transposing df1 and df2 by the column hour: i want to cluster usage per hour of the day.\n## transposing hour column\ndf1 = pd.pivot_table(df1, values = 'meter_reading', index=[\"building_id\", 'day'], columns = 'hour').reset_index()\ndf1 = df1.drop('day', 1)\n\n## filter rows where sum of 0 > 0 and the count of non null is > 23.\ndf1 = df1[df1.iloc[:,1:25].ne(0).sum(1) > 23 ]\n\n# filter rows where the sum of NaN from column 0 to 23 is less than 10\ndf1 = df1[df1.isnull().sum(axis=1) < 10] \n\n# fill NaN with row mean.\ndf1.iloc[:,1:25] = df1.iloc[:,1:25].T.fillna(df1.iloc[:,1:25].mean(axis=1)).T\n\n## compute the percentage of usage per hour per day.\n#df1.iloc[:,1:25] = 100 * df1.iloc[:,1:25].div(df1.iloc[:,1:25].sum(axis=1), axis=0) # 25776 * 25\n#df1 = df1.dropna(0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# df2\ndf2 = pd.pivot_table(df2, values = 'meter_reading', index=[\"building_id\", \"day\"], columns = 'hour').reset_index()\ndf2 = df2.drop('day', 1)\n\n## filter rows where sum of 0 > 0 and the count of non null is > 23.\ndf2 = df2[df2.iloc[:,1:25].ne(0).sum(1) > 23 ]\n\n# filter rows where the sum of NaN from column 0 to 23 is less than 10\ndf2 = df2[df2.isnull().sum(axis=1) < 10] \n\n# fill NaN with row mean.\ndf2.iloc[:,1:25] = df2.iloc[:,1:25].T.fillna(df2.iloc[:,1:25].mean(axis=1)).T\n\n## compute the percentage of usage per hour per day.\n#df2.loc[:,1:25] = 100 * df2.iloc[:,1:25].div(df2.iloc[:,1:25].sum(axis=1), axis=0) # 11349  * 25","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The script above does these:\n\n- pivot column hour: I am interested in visualizing the distribution of usage per hour while indentifying missing values.\n- remove all rows where there is no usage throughout the day as well as building where only one value per day is available; it is wrong data.\n- remove all rows where the count of missing values are more than 10. out of 24 hours, if 10 values are missing, remove rows.\n- filling remaining missing values with the average of the row per building.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#plot df1: weekdays\nplt.style.use('seaborn')\ndf1.iloc[:,1:25].T.plot(figsize=(16,8), legend=False, color='blue', alpha=0.01)\nplt.xlabel(\"hour of the day\")\nplt.ylabel(\"usage (in kwh)\")\nplt.title(\"Weekdays hourly usage\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.style.use('seaborn')\ndf2.iloc[:,1:25].T.plot(figsize=(16,8), legend=False, color='blue', alpha = 0.01)\nplt.xlabel(\"hour of the day\")\nplt.ylabel(\"usage (in kwh)\")\nplt.title(\"Weekends hourly usage\")\nplt.show()\n\n## possible clusters are grouped into categories of users (big: could be companies with running processes, small: households; regular consumers)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From both plots, we can observe:\n\n- weekdays: there seems to be 4 to 6 distinct groups: below 500, between 500 and more\n- weekends: there seems to be 4 to 6 groups (clusters): below 400, between 600 and 800, 1000 and 1400, above 1200\n\nWe can visually see this data can be clustered in type of users: low, medium, high and more. However, it is not possible to know for the optimum number of k by which data should be grouped and for which the error rate is the lowest. One way to find the optimum k is via the Elbow method which consists in:\n\n- Running kmeans clustering on different values ranging from 2 to 15 for example; for each value of k, the sum of squared errors (SSE) is calculated. \n- Plotting the line chart made of the number of potential k (x axis) and SSE (y axis). if the line ressembles like an arm, the optimum k is at the elbow.\n- The idea is the select small value of k that minimises SSE.\n ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Modelling: kmeans clustering.\n\nWhat is kmeans clustering and how does it work?\nIt is an unsupervised machine learning algorithm aiming at grouping data points based on similarities with the number of groups (clusters) represented by k. It follows these steps:\n\n- Choosing the best number of k.\n- Assigning each data point to the closest centroid by calculating its euclidian diantance with respoect to each centroid\n- Determine the new cluster centroid by computing the average of each cluster\n- Repeating steps 2 and 3 until none of the cluster assignemnts changes.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## weekdays\ndistortions = []\nK1 = range(1,10)\nfor k in K1:\n    kmeanModel1 = KMeans(n_clusters=k)\n    kmeanModel1.fit(df1.iloc[:, 1:25])\n    distortions.append(kmeanModel1.inertia_)\n\n# plot elbow line\nplt.figure(figsize=(16,8))\nplt.plot(K1, distortions)\nplt.xlabel('k')\nplt.ylabel('Distortion')\nplt.title('The Elbow Method showing the optimal k (weekdays)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The elbow line plot indicates optimum k is between 3 and 4. After 4, no improvement in SSE is observed. On this occasion, I will go for k = 4.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# weekdays model\nkmeans1 = KMeans(n_clusters=4).fit(df1.iloc[:, 1:25])\ncentroids1 = kmeans1.cluster_centers_\nprint(centroids1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## plot all weekdays centroids.\nplt.style.use('seaborn')\nax = sns.lineplot(x='hour', y='usage', marker=\"o\", \n                  hue=\"index\",palette=[\"C0\", \"C1\", \"C2\", \"C3\"],\n                  data=pd.melt(pd.DataFrame(centroids1).reset_index(), \n                               id_vars=\"index\", var_name=\"hour\", value_name=\"usage\")).set_title('clustering weekdays')\nplt.legend(title = 'Cluster', loc= 'upper right', labels = ['low', 'high', 'very high', 'medium'])  \nplt.show()                ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"- cluster low: can be household buildings where usage is very low\n- cluster medium: can be small businesses.\n- cluster high and very high: can be businesses with ongoing running processes at specific hours of the day. Both have similar patterns (at around 5am, usage starts to increase and peaks around 7am, gets constant until 15pm and then decreases.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"## weekends\ndistortions = []\nK2 = range(1,10)\nfor k in K2:\n    kmeanModel2 = KMeans(n_clusters=k)\n    kmeanModel2.fit(df2.iloc[:, 1:25])\n    distortions.append(kmeanModel2.inertia_)\n\n# plot elbow line\nplt.figure(figsize=(16,8))\nplt.plot(K2, distortions, 'bx-')\nplt.xlabel('k')\nplt.ylabel('Distortion')\nplt.title('The Elbow Method showing the optimal k (weekends)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# weekends model with k = 4\nkmeans2 = KMeans(n_clusters=4).fit(df2.iloc[:, 1:25])\ncentroids = kmeans2.cluster_centers_\nprint(centroids)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## plot all centroids.\nplt.style.use('seaborn')\nax = sns.lineplot(x='hour', y='usage', marker=\"o\", \n                  hue=\"index\",palette=[\"C0\", \"C1\", \"C2\", \"C3\"],\n                  data=pd.melt(pd.DataFrame(centroids).reset_index(), \n                               id_vars=\"index\", var_name=\"hour\", value_name=\"usage\")).set_title('clustering weekends')\nplt.legend(title = 'Cluster', loc= 'center right', labels = ['low', 'high', 'medium', 'very high'])  \nplt.show() ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"weekends clusters are similar to weekdays; however, the highest and very highest users groups comsume less electricity during weekends and the lines are flat as opposed to weekdays.\n\nFurther work: cluster each category; more pattern such as early morning, mid-morning, afternoon and night owl users  might be uncovered","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"","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":4}