{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Note\n\n- Please view in darkmode. \n\n- This kernel is inspired by Susan Li. Check her publications at https://towardsdatascience.com/@actsusanli"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd \n\nimport seaborn as sns \nimport matplotlib.dates as md\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.axes_grid1 import host_subplot\nimport mpl_toolkits.axisartist as AA\nplt.style.use(['fivethirtyeight', 'dark_background'])\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.cluster import KMeans\nfrom sklearn.covariance import EllipticEnvelope\nfrom sklearn.ensemble import IsolationForest\nfrom sklearn.svm import OneClassSVM\nfrom mpl_toolkits.mplot3d import Axes3D\n\n!pip install pyemma\nfrom pyemma import msm\n%matplotlib inline\n\nimport os\n\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/expedia-personalized-sort/data/train.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.isnull().sum()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Subset data\nSelect property / visitor location country / srch_room_count with the most data points "},{"metadata":{"trusted":true},"cell_type":"code","source":"# prop_id corresponding to \ntrain['prop_id'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train['visitor_location_country_id'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Num of rooms specified in search by customer\ntrain['srch_room_count'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Subset df \ndf = train.loc[train['prop_id'] == 104517]\n\ndf = df.loc[df['visitor_location_country_id'] == 219]\n\ndf = df.loc[df['srch_room_count'] == 1]\n\n# srch_saturday = if stay includes Sat night \n# srch_booking_window = num of days between search date and hotel stay start date \ndf = df[['date_time', 'price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Point anomaly = max price usd 5584"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.loc[(train['price_usd'] == 5584) & \n         (train['visitor_location_country_id'] == 219)]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Possible wrong search, no intention to book"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Remove 5584 \ndf = df.loc[df['price_usd'] < 5584]\ndf['price_usd'].describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['date_time'].min(), df['date_time'].max()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['date_time'].describe()\n\ndf['date_time'] = pd.to_datetime(df['date_time'])\n\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.plot(x = 'date_time', \n        y = 'price_usd', \n        figsize = (16, 8))\n\nplt.xlabel('dates')\nplt.ylabel('USD')\nplt.title('Time series of room price by date of search');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"a = df.loc[df['srch_saturday_night_bool'] == 0, 'price_usd']\nb = df.loc[df['srch_saturday_night_bool'] == 1, 'price_usd']\n\nplt.figure(figsize = (16, 8))\n\nplt.hist(a, bins = 80, \n         alpha = 0.3, \n         label = 'search w/o Sat night stay')\n\nplt.hist(b, bins = 80, \n         alpha = 0.3, \n         label = 'search w/ Sat night stay')\n\nplt.xlabel('Price')\nplt.ylabel('Freq')\nplt.legend()\nplt.title('Sat night search')\nplt.plot();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['srch_saturday_night_bool'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Kurtosis: %f' % df['price_usd'].kurt())\nprint('Skewness: %f' % df['price_usd'].skew())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(df['price_usd'], \n                 hist = False, label = 'USD')\n\nsns.distplot(df['srch_booking_window'], \n                  hist = False, label = 'booking window')\n\nplt.xlabel('dist')\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(a, hist = False, rug = False)\nsns.distplot(b, hist = False, rug = False)\n\nsns.despine()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df.sort_values('date_time')\ndf['date_time_int'] = df.date_time.astype(np.int64)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Cluster-based models\n\n- k-means\n- isolation forest\n- clustering\n\nPotential outliers:\n- usd\n- srch_booking_window (days between search and first stay date)\n- srch_saturday (stay includes sat night"},{"metadata":{},"cell_type":"markdown","source":"# K-means\n\n- create 'k' similar clusters of instances  \n- Instances outside of clusters = possible anomalies "},{"metadata":{"trusted":true},"cell_type":"code","source":"# Determine optimal cluster num using elbow method \ndata = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]\nn_cluster = range(1, 20)\n\nkmeans = [KMeans(n_clusters = i).fit(data) for i in n_cluster]\nscores = [kmeans[i].score(data) for i in range(len(kmeans))]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# elbow curve \nfig, ax = plt.subplots(figsize = (16, 8))\nax.plot(n_cluster, scores, color = 'orange')\n\nplt.xlabel('clusters num')\nplt.ylabel('score')\nplt.title('elbow curve')\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"set n_clusters to 7\n\nn_clusters > 7 = additional clusters do not explain greater variance in variable \n\nwhere variable = price_usd"},{"metadata":{"trusted":true},"cell_type":"code","source":"# k means output \nX = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]\nX = X.reset_index(drop = True)\n\nkm = KMeans(n_clusters = 7)\nkm.fit(X)\nkm.predict(X)\n\nlabels = km.labels_\n\nX.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3D clusters \nplot using k means output "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(1, figsize = (7, 7))\n\nax = Axes3D(fig, rect = [0, 0, 0.95, 1], \n            elev = 48, azim = 134)\n\nax.scatter(X.iloc[:, 0], \n           X.iloc[:, 1], \n           X.iloc[:, 2],\n           c = labels.astype(np.float), edgecolor = 'm')\n\nax.set_xlabel('USD')\nax.set_ylabel('srch_booking_window')\nax.set_zlabel('srch_saturday_night_bool')\n\nplt.title('K Means', fontsize = 10);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import pylab as pl ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PCA "},{"metadata":{"trusted":true},"cell_type":"code","source":"#Y = df[['price_usd']]\n#X = df[['srch_booking_window']]\n\n#Nc = range(1, 20)\n#kmeans = [KMeans(n_clusters = i) for i in Nc]\n\n#score = [kmeans[i].fit(Y).score(Y) for i in range(len(kmeans))]\n\n#plt.figure(figsize = (16, 8))\n#pl.plot(Nc, score)\n#pl.xlabel('cluster num')\n#pl.ylabel('score')\n#pl.title('elbow curve')\n#pl.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#pca = PCA(n_components = 1).fit(Y)\n\n#pca_d = pca.transform(Y)\n#pca_c = pca.transform(X)\n\n#kmeans = KMeans(n_clusters = 7)\n#kmeansoutput = kmeans.fit(Y)\n\n#pl.figure('7 cluster k-means')\n#pl.figure(figsize = (16, 8))\n\n#pl.scatter(pca_c[:, 0], \n#           pca_d[:, 0], \n#           c = kmeansoutput.labels_)\n\n#pl.xlabel('booking window')\n#pl.ylabel('USD')\n#pl.title('7 cluster')\n#pl.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]\n\nX = data.values\nX_std = StandardScaler().fit_transform(X)\n\n# Calc eigenvec cor & eig_vals of covar matrix \nmean_vec = np.mean(X_std, axis = 0)\n\ncov_mat = np.cov(X_std.T)\n\neig_vals, eig_vecs = np.linalg.eig(cov_mat)\n\n# eig_val,eig_vecs tuple\neig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]\n\neig_pairs.sort(key = lambda x: x[0], reverse = True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Calc explained var from eig_vals \ntotal = sum(eig_vals)\n\n# Individual explained var \nvar_exp = [(i/total)*100 for i in sorted(eig_vals, reverse = True)]\n\n# Cumulative explained var \ncum_var_exp = np.cumsum(var_exp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (16, 8))\nplt.bar(range(len(var_exp)), var_exp, \n        alpha = 0.5, align = 'center', \n        label = 'individual explained var', \n        color = 'r'\n       )\n\nplt.step(range(len(cum_var_exp)), cum_var_exp,\n         where = 'mid',\n         label = 'cumulative explained var')\n\nplt.xlabel('principal components')\nplt.ylabel('explained var ratio')\nplt.legend(loc = 'best')\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Component 1 explains approx 50% of var. \n2 = explains < 40\n3 = explains < 20 \n\nComponents 1 + 2 = explain approx 80% of var \n\n- set n_components = 2\n- standardize features "},{"metadata":{"trusted":true},"cell_type":"code","source":"data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]\n\n# Standardize features\nX_std = StandardScaler().fit_transform(X)\ndata = pd.DataFrame(X_std)\n\n# Reduce components to 2 \npca = PCA(n_components = 2)\ndata = pca.fit_transform(data)\n\n# Standardize 2 new features \nscaler = StandardScaler()\nnp_scaled = scaler.fit_transform(data)\ndata = pd.DataFrame(np_scaled)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"kmeans = [KMeans(n_clusters = i).fit(data) for i in n_cluster]\n\ndf['cluster'] = kmeans[7].predict(data)\ndf.index = data.index\n\ndf['pc1'] = data[0]\ndf['pc2'] = data[1]\ndf['cluster'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def getDistanceByPoint(data, model):\n    distance = pd.Series()\n    for i in range(0,len(data)):\n        Xa = np.array(data.loc[i])\n        Xb = model.cluster_centers_[model.labels_[i]-1]\n        distance.set_value(i, np.linalg.norm(Xa-Xb))\n    return distance","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"outliers_fraction = 0.01\n\ndistance = getDistanceByPoint(data, kmeans[9])\noutlier_num = int(outliers_fraction * len(distance))\n\nthreshold = distance.nlargest(outlier_num).min()\n\ndf['anomaly'] = (distance >= threshold).astype(int)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (12, 6))\n\ncolors = {0:'blue', 1:'red'}\n\nax.scatter(df['pc1'], df['pc2'], \n           c = df['anomaly'].apply(lambda x: colors[x]))\n\nplt.xlabel('pc1')\nplt.ylabel('pc2')\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df.sort_values('date_time')\ndf['date_time'] = df.date_time.astype(np.int64)\n\n# object with anomalies\na = df.loc[df['anomaly'] == 1, \n           ['date_time_int', 'price_usd']]\n\na","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (10, 5))\n\nax.plot(df['date_time_int'], df['price_usd'], \n        color = 'orange', label = 'Normal')\n\nax.scatter(a['date_time_int'], a['price_usd'],\n           color = 'red', label = 'Anomaly')\n\nplt.xlabel('time')\nplt.ylabel('USD')\nplt.legend()\nplt.show();\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['anomaly'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"a = df.loc[df['anomaly'] == 0, 'price_usd']\nb = df.loc[df['anomaly'] == 1, 'price_usd']\n\nfig, axs = plt.subplots(figsize = (10, 5))\naxs.hist([a, b], \n         bins = 50, stacked = True, \n         color = ['orange', 'red'])\n\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.anomaly.value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Isolation Forest"},{"metadata":{},"cell_type":"markdown","source":"Detect anomalies based on data points that are few and different \n\n- No use of density / distance measure \n    i.e. different from clustering based / distanced based algorithms \n\n- Randomly select a feature \n    \n- Randomly select a split between max and min values of selected feature \n    \n- Length of path, avged over a forest of random trees = measure of normality \n\n- Random partitioning = shorter path for anomalies\n\n- If forest produces shorter paths for samples, then they are likely to be anomalies "},{"metadata":{"trusted":true},"cell_type":"code","source":"data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]\n\nscaler = StandardScaler()\nnp_scaled = scaler.fit_transform(data)\n\ndata = pd.DataFrame(np_scaled)\n\n# Isolation forest \noutliers_fraction = 0.01\nifo = IsolationForest(contamination = outliers_fraction)\n\nifo.fit(data)\n\ndf['anomaly1'] = pd.Series(ifo.predict(data))\n\nfig, ax = plt.subplots(figsize = (10, 5))\n\na = df.loc[df['anomaly1'] == -1, ['date_time_int', 'price_usd']]\n\nax.plot(df['date_time_int'], df['price_usd'], \n        color = 'orange', label = 'Normal')\n\nax.scatter(a['date_time_int'], a['price_usd'], \n           color = 'red', label = 'Anomaly')\n\nplt.legend()\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['anomaly1'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"a = df.loc[df['anomaly1'] == 1, 'price_usd']\nb = df.loc[df['anomaly1'] == -1, 'price_usd']\n\nfig, ax = plt.subplots(figsize = (10, 5))\n\nax.hist([a, b],\n        bins = 50, stacked = True, \n        color = ['orange', 'red'] )\n\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['anomaly1'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Support vector machine models\n\nAssociated with supervised learning \n\n- One class SVM\n- Gaussian dist \n- Markov chain"},{"metadata":{},"cell_type":"markdown","source":"# One class SVM \n- unsupervised anomaly detection \n- estimate support of high dimensional distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]\nscaler = StandardScaler()\nnp_scaled = scaler.fit_transform(data)\n\ndata = pd.DataFrame(np_scaled)\n\n# Train \n\nosvm = OneClassSVM(nu = outliers_fraction, \n                   kernel = 'rbf', \n                   gamma = 0.01)\n\nosvm.fit(data)\n\ndf['anomaly2'] = pd.Series(osvm.predict(data))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize = (10, 5))\n\na = df.loc[df['anomaly2'] == -1, \n           ['date_time_int', 'price_usd']]\n\nax.plot(df['date_time_int'], df['price_usd'], \n        color = 'orange', \n        label = 'Normal')\n\nax.scatter(a['date_time_int'], a['price_usd'], \n           color = 'red', \n           label = 'Anomaly')\n\nplt.legend()\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"a = df.loc[df['anomaly2'] == 1, 'price_usd']\nb = df.loc[df['anomaly2'] == -1, 'price_usd']\n\nfig, ax = plt.subplots(figsize = (10, 5))\n\nax.hist([a, b], bins = 50, \n        stacked = True, color = ['orange','red'])\n\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['anomaly2'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Gaussian distribution\n\n- Assume data is normally distributed\n\n- Use covariance.EllipticEnvelope from scikit-learn to find key params of general distribution by assuming entire dataset = an expression of an underlying multivariate Gaussian distribution\n"},{"metadata":{},"cell_type":"markdown","source":"Create two dfs based on categories defined by sat boolean"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_class0 = df.loc[df['srch_saturday_night_bool'] == 0, 'price_usd']\ndf_class1 = df.loc[df['srch_saturday_night_bool'] == 1, 'price_usd']\n\nfig, axs = plt.subplots(1, 2)\n\ndf_class0.hist(ax = axs[0], bins = 50, color = 'orange')\ndf_class1.hist(ax = axs[1], bins = 50, color = 'red');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Apply EllipticEnvelope to each category \n\nSet contamination param (proportion of outliers present in dataset)\n\nUse decision function to compute decision function of given observations (equivalent to shifted Mahalanobis distances. \n\nThreshold for identifying as outliers = 0 (compatible with other detection algorithms)\n\npredict(x_train) predict labels of X_train according to fitted model\n\n1 = normal\n-1 = anomaly"},{"metadata":{"trusted":true},"cell_type":"code","source":"envelope = EllipticEnvelope(contamination = outliers_fraction)\n\nx_train = df_class0.values.reshape(-1, 1)\nenvelope.fit(x_train)\n\ndf_class0 = pd.DataFrame(df_class0)\ndf_class0['deviation'] = envelope.decision_function(x_train)\ndf_class0['anomaly'] = envelope.predict(x_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"envelope = EllipticEnvelope(contamination = outliers_fraction)\n\nx_train = df_class1.values.reshape(-1, 1)\nenvelope.fit(x_train)\n\ndf_class1 = pd.DataFrame(df_class1)\ndf_class1['deviation'] = envelope.decision_function(x_train)\ndf_class1['anomaly'] = envelope.predict(x_train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_class = pd.concat([df_class0, df_class1])\ndf['anomaly3'] = df_class['anomaly']\n\nfig, ax = plt.subplots(figsize = (10, 5))\n\na = df.loc[df['anomaly3'] == -1, \n           ('date_time_int', 'price_usd')]\n\nax.plot(df['date_time_int'], df['price_usd'], \n        color = 'orange')\n\nax.scatter(a['date_time_int'], a['price_usd'],\n          color = 'red')\n\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['anomaly3'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"a = df.loc[df['anomaly3'] == 1, 'price_usd']\nb = df.loc[df['anomaly3'] == -1, 'price_usd']\n\nfig, ax = plt.subplots(figsize = (10, 5))\nax.hist([a, b], \n        bins = 50, stacked = True, \n        color = ['orange', 'red'])\n\nplt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Results \nanomalies detected only show abnormally high prices, no abnormally low prices "}],"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}