{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# New York City Taxi Trip Duration","metadata":{"_uuid":"8b803d21f9fc2862a5f53d730c80f8f045fe1b7c"}},{"cell_type":"markdown","source":"The competition dataset is based on the 2016 NYC Yellow Cab trip record data made available in Big Query on Google Cloud Platform. The data was originally published by the NYC Taxi and Limousine Commission (TLC). The data was sampled and cleaned for the purposes of this playground competition. Based on individual trip attributes, participants should **predict the duration of each trip in the test set.**","metadata":{"_uuid":"928ed469171a0ae01e24043f2753b488b59f205f","_kg_hide-input":false}},{"cell_type":"markdown","source":"## I. Data loading and overview <a id=\"one\"></a>","metadata":{"_uuid":"8bcf7a799ff0f92329de70fff29e31c9a62d97cd"}},{"cell_type":"code","source":"import numpy as np\nnp.random.seed(42)\nimport pandas as pd\n\nimport shap\n\nfrom sklearn.metrics import mean_squared_log_error as MSE\nfrom sklearn.model_selection import train_test_split\n\nimport xgboost as xgb\nimport lightgbm as lgb\n","metadata":{"_uuid":"a2d92f36f38765c5fbed4d672ea8f746cb073202","execution":{"iopub.status.busy":"2022-03-29T14:40:35.335476Z","iopub.execute_input":"2022-03-29T14:40:35.33571Z","iopub.status.idle":"2022-03-29T14:40:39.782206Z","shell.execute_reply.started":"2022-03-29T14:40:35.335647Z","shell.execute_reply":"2022-03-29T14:40:39.781492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### a. Loading the data <a id=\"one-a\"></a>","metadata":{"_uuid":"98d81f72eded731792db8202cd9ffb7fb2a4c6c2"}},{"cell_type":"code","source":"df = pd.read_csv(\"/kaggle/input/nyc-taxi-trip-duration/train.zip\")\ntest = pd.read_csv(\"/kaggle/input/nyc-taxi-trip-duration/test.zip\")","metadata":{"_uuid":"573a8353e37dfbfd0c193a7118199528c38650f7","execution":{"iopub.status.busy":"2022-03-29T14:40:39.827304Z","iopub.execute_input":"2022-03-29T14:40:39.827568Z","iopub.status.idle":"2022-03-29T14:40:46.837871Z","shell.execute_reply.started":"2022-03-29T14:40:39.827535Z","shell.execute_reply":"2022-03-29T14:40:46.836925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"add_on_df = pd.read_csv(\"/kaggle/input/2014-new-york-city-taxi-trips/nyc_taxi_data_2014.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:46.853944Z","iopub.execute_input":"2022-03-29T14:40:46.854374Z","iopub.status.idle":"2022-03-29T14:40:46.858607Z","shell.execute_reply.started":"2022-03-29T14:40:46.854337Z","shell.execute_reply":"2022-03-29T14:40:46.857729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"add_on_df.dropna(inplace=True,)\ntrip_duration = pd.to_datetime(add_on_df.dropoff_datetime) - pd.to_datetime(add_on_df.pickup_datetime)\ntrip_duration /= np.timedelta64(1, 's')\nadd_on_df[\"trip_duration\"] = trip_duration\nadd_on_df.drop(add_on_df[add_on_df.trip_duration < 0].index, inplace=True)\ntrip_duration = add_on_df[\"trip_duration\"]\n\nrevelant_columns = list(test.columns) \nrevelant_columns.remove(\"id\")\n\nadd_on_df = add_on_df[revelant_columns]\n# add_on_df[\"vendor_id\"] = add_on_df[\"vendor_id\"].astype('category')\n# add_on_df[\"vendor_id\"] = add_on_df[\"vendor_id\"].cat.codes\nadd_on_df.shape, trip_duration.shape","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:46.859955Z","iopub.execute_input":"2022-03-29T14:40:46.860698Z","iopub.status.idle":"2022-03-29T14:40:46.867205Z","shell.execute_reply.started":"2022-03-29T14:40:46.860662Z","shell.execute_reply":"2022-03-29T14:40:46.866472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# To do Generalisation Hypothesis Experiments, uncomment the following to assign the add-on df to test variable, overwriting the testing test\n# test = add_on_df[:250_000]\n# trip_duration = trip_duration[:250_000]","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:46.880606Z","iopub.execute_input":"2022-03-29T14:40:46.88091Z","iopub.status.idle":"2022-03-29T14:40:46.887927Z","shell.execute_reply.started":"2022-03-29T14:40:46.880876Z","shell.execute_reply":"2022-03-29T14:40:46.887177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test.shape","metadata":{"_uuid":"bafeb9384f26ee7c8f3ac08aea33991636490b61","execution":{"iopub.status.busy":"2022-03-29T14:40:46.892826Z","iopub.execute_input":"2022-03-29T14:40:46.893028Z","iopub.status.idle":"2022-03-29T14:40:46.898561Z","shell.execute_reply.started":"2022-03-29T14:40:46.893004Z","shell.execute_reply":"2022-03-29T14:40:46.897767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### b. Overview <a id=\"one-b\"></a>","metadata":{"_uuid":"1d42ab54483a991704abe137e8ada2ec839484a7"}},{"cell_type":"markdown","source":"Column | Description\n------- | -------\n**id** | a unique identifier for each trip  \n**vendor_id** | a code indicating the provider associated with the trip record  \n**pickup_datetime** | date and time when the meter was engaged  \n**dropoff_datetime** | date and time when the meter was disengaged  \n**passenger_count** | the number of passengers in the vehicle (driver entered value)  \n**pickup_longitude** | the longitude where the meter was engaged  \n**pickup_latitude** | the latitude where the meter was engaged  \n**dropoff_longitude** | the longitude where the meter was disengaged  \n**dropoff_latitude** | the latitude where the meter was disengaged  \n**store_and_fwd_flag** | This flag indicates whether the trip record was held in vehicle memory before sending to the vendor because the vehicle did not have a connection to the server (Y=store and forward; N=not a store and forward trip)  \n**trip_duration** | duration of the trip in seconds  \n\n*Disclaimer: The decision was made to not remove dropoff coordinates from the dataset order to provide an expanded set of variables to use in Kernels.*","metadata":{"_uuid":"cc9386ab0a11c674223c6ecc3604085837f7ce77"}},{"cell_type":"markdown","source":"## II. Data cleaning <a id=\"two\"></a>","metadata":{"_uuid":"d85baac0eff4e6525df407b262596a72e8a5fecb"}},{"cell_type":"markdown","source":"### a. Duplicated and missing values <a id=\"two-a\"></a>","metadata":{"_uuid":"d74e5da2bb632a78ada91ef797bf1a562e9df80c"}},{"cell_type":"code","source":"#Count the number of duplicated rows\ndf.duplicated().sum()","metadata":{"_uuid":"b4da4373e43625132d7df1723ed12c8253b2a007","execution":{"iopub.status.busy":"2022-03-29T14:40:47.792015Z","iopub.execute_input":"2022-03-29T14:40:47.792842Z","iopub.status.idle":"2022-03-29T14:40:49.756212Z","shell.execute_reply.started":"2022-03-29T14:40:47.792803Z","shell.execute_reply":"2022-03-29T14:40:49.755534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Count the number of NaN values for each column\ndf.isna().sum()","metadata":{"_uuid":"daf4dc10c14b5191a0a5b9d4576bcab57160f726","execution":{"iopub.status.busy":"2022-03-29T14:40:49.757526Z","iopub.execute_input":"2022-03-29T14:40:49.757783Z","iopub.status.idle":"2022-03-29T14:40:50.322931Z","shell.execute_reply.started":"2022-03-29T14:40:49.757747Z","shell.execute_reply":"2022-03-29T14:40:50.322223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# add_on_df.isna().sum()","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:50.324394Z","iopub.execute_input":"2022-03-29T14:40:50.324678Z","iopub.status.idle":"2022-03-29T14:40:50.328178Z","shell.execute_reply.started":"2022-03-29T14:40:50.324643Z","shell.execute_reply":"2022-03-29T14:40:50.327361Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# add_on_df.shape","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:50.330607Z","iopub.execute_input":"2022-03-29T14:40:50.33155Z","iopub.status.idle":"2022-03-29T14:40:50.338484Z","shell.execute_reply.started":"2022-03-29T14:40:50.33151Z","shell.execute_reply":"2022-03-29T14:40:50.337814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are no duplicated or missing values.","metadata":{"_uuid":"d73b386029446bcde76f881fdbb5357d4f08964c"}},{"cell_type":"markdown","source":"### b. Deal with outliers <a id=\"two-b\"></a>","metadata":{"_uuid":"810ae616deaec1e0b1ea42e7bf90a253e6676ce4"}},{"cell_type":"markdown","source":"We clearly see `trip_duration` takes strange values for `min` and `max`. Let's have a quick visualization with a boxplot.","metadata":{"_uuid":"72e367455010f3ca0ad92a478888ad9a6c8d7fcc"}},{"cell_type":"markdown","source":"There are outliers for `trip_duration`. I can't find a proper interpretation and it will probably damage our model, so I choose to get rid of them.","metadata":{"_uuid":"1b7d073d9b3954eb84c1d3ca0337c0f488c2c88a"}},{"cell_type":"code","source":"#Only keep trips that lasted less than 5900 seconds\ndf = df[(df.trip_duration < 5900)]\n\n#Only keep trips with passengers\ndf = df[(df.passenger_count > 0)]","metadata":{"_uuid":"716e1c857c7e526009cb09217cf9a929ded2bd34","execution":{"iopub.status.busy":"2022-03-29T14:40:50.67219Z","iopub.execute_input":"2022-03-29T14:40:50.672633Z","iopub.status.idle":"2022-03-29T14:40:50.830124Z","shell.execute_reply.started":"2022-03-29T14:40:50.672593Z","shell.execute_reply":"2022-03-29T14:40:50.829383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## III. Features engineering <a id=\"three\"></a>","metadata":{"_uuid":"ded8fe370f88295307313d5ccb1d6c909e573f59"}},{"cell_type":"markdown","source":"### b. Deal with categorical features <a id=\"three-b\"></a>","metadata":{"_uuid":"46ba9f740e659e6cca0a749e0a2288d48e3fc78d"}},{"cell_type":"code","source":"#One-hot encoding binary categorical features\ndf = pd.concat([df, pd.get_dummies(df['store_and_fwd_flag'])], axis=1)\ntest = pd.concat([test, pd.get_dummies(test['store_and_fwd_flag'])], axis=1)\n\ndf.drop(['store_and_fwd_flag'], axis=1, inplace=True)\ntest.drop(['store_and_fwd_flag'], axis=1, inplace=True)\n\n\ndf = pd.concat([df, pd.get_dummies(df['vendor_id'])], axis=1)\ntest = pd.concat([test, pd.get_dummies(test['vendor_id'])], axis=1)\n\ndf.drop(['vendor_id'], axis=1, inplace=True)\ntest.drop(['vendor_id'], axis=1, inplace=True)","metadata":{"_uuid":"84f46b633163790904b2c1f9476943926d64b272","execution":{"iopub.status.busy":"2022-03-29T14:40:50.996303Z","iopub.execute_input":"2022-03-29T14:40:50.996564Z","iopub.status.idle":"2022-03-29T14:40:52.421031Z","shell.execute_reply.started":"2022-03-29T14:40:50.99653Z","shell.execute_reply":"2022-03-29T14:40:52.420296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"ym <a id=\"three-c\"></a>","metadata":{"_uuid":"727031c83cb7bbd8b6f30ccfb9b72f12dfebcfd2"}},{"cell_type":"code","source":"#Datetyping the dates\ndf['pickup_datetime'] = pd.to_datetime(df.pickup_datetime)\ntest['pickup_datetime'] = pd.to_datetime(test.pickup_datetime)\n\ndf.drop(['dropoff_datetime'], axis=1, inplace=True) #as we don't have this feature in the testset\n\n#Date features creations and deletions\ndf['month'] = df.pickup_datetime.dt.month\ndf['week'] = df.pickup_datetime.dt.week\ndf['weekday'] = df.pickup_datetime.dt.weekday\ndf['hour'] = df.pickup_datetime.dt.hour\ndf['minute'] = df.pickup_datetime.dt.minute\ndf['minute_oftheday'] = df['hour'] * 60 + df['minute']\ndf.drop(['minute'], axis=1, inplace=True)\n\ntest['month'] = test.pickup_datetime.dt.month\ntest['week'] = test.pickup_datetime.dt.week\ntest['weekday'] = test.pickup_datetime.dt.weekday\ntest['hour'] = test.pickup_datetime.dt.hour\ntest['minute'] = test.pickup_datetime.dt.minute\ntest['minute_oftheday'] = test['hour'] * 60 + test['minute']\ntest.drop(['minute'], axis=1, inplace=True)\n\ndf.drop(['pickup_datetime'], axis=1, inplace=True)\n\ndf.info()","metadata":{"scrolled":true,"_uuid":"7b625def5299b397333dd0c900cfdd50543bc6c8","execution":{"iopub.status.busy":"2022-03-29T14:40:52.422464Z","iopub.execute_input":"2022-03-29T14:40:52.422716Z","iopub.status.idle":"2022-03-29T14:40:55.063112Z","shell.execute_reply.started":"2022-03-29T14:40:52.422682Z","shell.execute_reply":"2022-03-29T14:40:55.062355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### d. Distance and speed creations <a id=\"three-d\"></a>","metadata":{"_uuid":"82a18eb9e511ace5bc2eedeff3fac70e0860b0a0"}},{"cell_type":"code","source":"#Function aiming at calculating distances from coordinates\ndef ft_haversine_distance(lat1, lng1, lat2, lng2):\n    lat1, lng1, lat2, lng2 = map(np.radians, (lat1, lng1, lat2, lng2))\n    AVG_EARTH_RADIUS = 6371 #km\n    lat = lat2 - lat1\n    lng = lng2 - lng1\n    d = np.sin(lat * 0.5) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(lng * 0.5) ** 2\n    h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d))\n    return h\n\n#Add distance feature\ndf['distance'] = ft_haversine_distance(df['pickup_latitude'].values,\n                                                 df['pickup_longitude'].values, \n                                                 df['dropoff_latitude'].values,\n                                                 df['dropoff_longitude'].values)\ntest['distance'] = ft_haversine_distance(test['pickup_latitude'].values, \n                                                test['pickup_longitude'].values, \n                                                test['dropoff_latitude'].values, \n                                                test['dropoff_longitude'].values)","metadata":{"_uuid":"6e147135c6762619166754a4a16f380210513f61","execution":{"iopub.status.busy":"2022-03-29T14:40:55.064273Z","iopub.execute_input":"2022-03-29T14:40:55.065586Z","iopub.status.idle":"2022-03-29T14:40:55.212229Z","shell.execute_reply.started":"2022-03-29T14:40:55.065542Z","shell.execute_reply":"2022-03-29T14:40:55.21148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Function aiming at calculating the direction\ndef ft_degree(lat1, lng1, lat2, lng2):\n    AVG_EARTH_RADIUS = 6371 #km\n    lng_delta_rad = np.radians(lng2 - lng1)\n    lat1, lng1, lat2, lng2 = map(np.radians, (lat1, lng1, lat2, lng2))\n    y = np.sin(lng_delta_rad) * np.cos(lat2)\n    x = np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(lng_delta_rad)\n    return np.degrees(np.arctan2(y, x))\n\n#Add direction feature\ndf['direction'] = ft_degree(df['pickup_latitude'].values,\n                                df['pickup_longitude'].values,\n                                df['dropoff_latitude'].values,\n                                df['dropoff_longitude'].values)\ntest['direction'] = ft_degree(test['pickup_latitude'].values,\n                                  test['pickup_longitude'].values, \n                                  test['dropoff_latitude'].values,\n                                  test['dropoff_longitude'].values)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:55.21365Z","iopub.execute_input":"2022-03-29T14:40:55.21391Z","iopub.status.idle":"2022-03-29T14:40:55.476694Z","shell.execute_reply.started":"2022-03-29T14:40:55.213874Z","shell.execute_reply":"2022-03-29T14:40:55.475937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Remove distance outliers\ndf = df[(df.distance < 200)]\n\n#Create speed feature\ndf['speed'] = df.distance / df.trip_duration\n\n#Remove speed outliers\ndf = df[(df.speed < 30)]\ndf.drop(['speed'], axis=1, inplace=True)","metadata":{"_uuid":"f085828ddbe1b1cd4e8dc1ab49874c3f43399810","execution":{"iopub.status.busy":"2022-03-29T14:40:55.478128Z","iopub.execute_input":"2022-03-29T14:40:55.478386Z","iopub.status.idle":"2022-03-29T14:40:55.654454Z","shell.execute_reply.started":"2022-03-29T14:40:55.478351Z","shell.execute_reply":"2022-03-29T14:40:55.653691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## IV. Model selection <a id=\"four\"></a>","metadata":{"_uuid":"b51df1f4338fa4128cbab12d8459bcf0c42728e6"}},{"cell_type":"markdown","source":"### a. Split <a id=\"four-a\"></a>\n\nthe 250_000 train_size 250_000 test size","metadata":{"_uuid":"5518276bcaa6c2ba0b13bf9fc3844e7b031235f9"}},{"cell_type":"code","source":"#Split the labeled data frame into two sets: features and target\ny = df[\"trip_duration\"]\ndf.drop([\"trip_duration\"], axis=1, inplace=True)\ndf.drop(['id'], axis=1, inplace=True)\nX = df\n\nX.shape, y.shape","metadata":{"_uuid":"4618aab947bb372d1c367ec7c12235f1f4f446c4","execution":{"iopub.status.busy":"2022-03-29T14:40:55.93963Z","iopub.execute_input":"2022-03-29T14:40:55.939903Z","iopub.status.idle":"2022-03-29T14:40:56.109991Z","shell.execute_reply.started":"2022-03-29T14:40:55.939868Z","shell.execute_reply":"2022-03-29T14:40:56.109311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Split the labeled data frame into two sets to train then test the models\n\nX_train, X_rest, y_train, y_rest = train_test_split(X, y, train_size=250_000, random_state=42)\n\nX_test, X_rest, y_test, y_rest = train_test_split(X_rest, y_rest, train_size=250_000, random_state=42)\n\nX_train.shape, y_train.shape, ","metadata":{"_uuid":"cfa6f955aafef42aaccc26674694cba2defa1680","execution":{"iopub.status.busy":"2022-03-29T14:40:56.126565Z","iopub.execute_input":"2022-03-29T14:40:56.127061Z","iopub.status.idle":"2022-03-29T14:40:56.725018Z","shell.execute_reply.started":"2022-03-29T14:40:56.127006Z","shell.execute_reply":"2022-03-29T14:40:56.724298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_on_df_clean_up():\n    test.drop('pickup_datetime', axis=1, inplace=True)\n    test.insert(8,2,0)\n    test.rename(columns={\"CMT\": 1},inplace=True)\n    X_add_on, _, y_add_on, _ = train_test_split(test, trip_duration, train_size=75_000, random_state=42)\n    return X_add_on, y_add_on","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Uncomment the following for Generalisation Hypothesis tests\n# X_add_on, y_add_on = add_on_df_clean_up()","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:56.787247Z","iopub.execute_input":"2022-03-29T14:40:56.788144Z","iopub.status.idle":"2022-03-29T14:40:56.793519Z","shell.execute_reply.started":"2022-03-29T14:40:56.788108Z","shell.execute_reply":"2022-03-29T14:40:56.792787Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### b. Metrics <a id=\"four-b\"></a>\nFor this specific problem, we'll measure the error using the RMSLE (Root Mean Squared Log Error).","metadata":{"_uuid":"dd61fbcbadd97bf9f9ec4fb198fb1138a2e8aaee"}},{"cell_type":"markdown","source":"### c. Models <a id=\"four-c\"></a>","metadata":{"_uuid":"c1b4e15a839bbf7eb394987c122956a8e3192b33"}},{"cell_type":"markdown","source":"#### classic gradient boosting\n\n0.7489497458840355 0.7434332056656134  \nCPU times: user 20.9 s, sys: 76.7 ms, total: 21 s  \nWall time: 21 s  \n0.4355758648932763  \ndoesn't perfmorm well enough to be used","metadata":{}},{"cell_type":"code","source":"# %%time\n# #Try GradientBoosting\n# from sklearn.ensemble import GradientBoostingRegressor\n\n# gb = GradientBoostingRegressor()\n# gb.fit(X_train, y_train)\n# print(gb.score(X_train, y_train), gb.score(X_test, y_test))\n\n# print(np.sqrt(MSE(y_test, np.clip(gb.predict(X_test),0,None))))","metadata":{"_uuid":"683114022fd7afcac599767973940ea30a5261d4","execution":{"iopub.status.busy":"2022-03-29T14:40:56.795279Z","iopub.execute_input":"2022-03-29T14:40:56.796267Z","iopub.status.idle":"2022-03-29T14:40:56.80232Z","shell.execute_reply.started":"2022-03-29T14:40:56.796224Z","shell.execute_reply":"2022-03-29T14:40:56.801586Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### random forest \n**with full data 80/20 split**  \n0.876998567217304 0.813376227423164  \nCPU times: user 7min 28s, sys: 609 ms, total: 7min 29s  \nWall time: 7min 31s  \n0.3764782338436293\n","metadata":{}},{"cell_type":"code","source":"# %%time\n# #Try RandomForest\n# from sklearn.ensemble import RandomForestRegressor\n\n# rfm = RandomForestRegressor(bootstrap=True,max_depth=90,max_features='auto',min_samples_split=15,min_samples_leaf=10,n_estimators=30)\n# rfm.fit(X_train, y_train)\n# # print(rfm.score(X_train, y_train), rfm.score(X_test, y_test),rfm.score(test, trip_duration))\n","metadata":{"_uuid":"79a3c5323aaa5259f2d170057a62e29c14582dd4","execution":{"iopub.status.busy":"2022-03-29T14:40:56.803415Z","iopub.execute_input":"2022-03-29T14:40:56.804128Z","iopub.status.idle":"2022-03-29T14:40:56.810651Z","shell.execute_reply.started":"2022-03-29T14:40:56.804089Z","shell.execute_reply":"2022-03-29T14:40:56.809938Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Light GBM Models\n**with full data 80/20 split**  \n\n\n**model 1:**  \nCPU times: user 1min 12s, sys: 1.07 s, total: 1min 13s  \nWall time: 39.8 s  \n0.3647868165914379  \n\n\n**model 2:**  \n0.9220040551405524 0.8444362948661079  \nCPU times: user 5min 26s, sys: 3.87 s, total: 5min 30s  \nWall time: 2min 55s  \n0.36742497049051626  ","metadata":{}},{"cell_type":"code","source":"%%time\n\nlgb_params = {\n    'metric': 'rmse',\n    'is_training_metric': True,\n    'learning_rate': 0.1,\n    'max_depth': 25,\n    'num_leaves': 1000, \n    'objective': 'regression',\n    'feature_fraction': 0.9,\n    'bagging_fraction': 0.5,\n    'max_bin': 1000 ,\n}\n\nlgb_train = lgb.Dataset(X_train, y_train)\nlgb_test = lgb.Dataset(X_test, y_test)\nlgbm1 = lgb.train(lgb_params, lgb_train, num_boost_round=100, valid_sets=[lgb_train], early_stopping_rounds=5,verbose_eval = -1)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:40:56.812709Z","iopub.execute_input":"2022-03-29T14:40:56.813538Z","iopub.status.idle":"2022-03-29T14:41:10.388506Z","shell.execute_reply.started":"2022-03-29T14:40:56.8135Z","shell.execute_reply":"2022-03-29T14:41:10.387932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# print(lgbm1.score(X_train, y_train), lgbm1.score(X_test, y_test), lgbm1.score(test, trip_duration))","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:41:10.390922Z","iopub.execute_input":"2022-03-29T14:41:10.391265Z","iopub.status.idle":"2022-03-29T14:41:10.395226Z","shell.execute_reply.started":"2022-03-29T14:41:10.391228Z","shell.execute_reply":"2022-03-29T14:41:10.394597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nfrom lightgbm import LGBMRegressor\n\nlgbm2 = lgb.LGBMRegressor(n_estimators=500, num_leaves=1000, max_depth=25, objective='regression') #0.39507740906150585\n# lgbm = lgb.LGBMRegressor() #0.3990924865808348\nlgbm2.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:41:10.397306Z","iopub.execute_input":"2022-03-29T14:41:10.397564Z","iopub.status.idle":"2022-03-29T14:41:41.881419Z","shell.execute_reply.started":"2022-03-29T14:41:10.39753Z","shell.execute_reply":"2022-03-29T14:41:41.880658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# print(lgbm2.score(X_train, y_train), lgbm2.score(X_test, y_test), lgbm2.score(test, trip_duration))\n# print(lgbm2.score(test, trip_duration[:250_000]))\n\n#Output:\n    #0.7812886118508641 0.7827256176145024\n    #0.3623481127815768\n    #CPU times: user 42 s, sys: 1.08 s, total: 43 s\n    #Wall time: 22.5 s","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:41:41.882616Z","iopub.execute_input":"2022-03-29T14:41:41.882937Z","iopub.status.idle":"2022-03-29T14:41:41.88878Z","shell.execute_reply.started":"2022-03-29T14:41:41.882899Z","shell.execute_reply":"2022-03-29T14:41:41.888081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### XGB\n\n**with full data 80/20 split**  \n0.9088948471799362 0.8339765297679975  \nCPU times: user 9min 33s, sys: 1.89 s, total: 9min 35s  \nWall time: 5min 25s   \n0.3674661170369321","metadata":{}},{"cell_type":"code","source":"xgb_params = {\n    'booster':            'gbtree',\n    'objective':          'reg:squarederror',\n    'learning_rate':      0.1,\n    'max_depth':          14,\n    'subsample':          0.8,\n    'colsample_bytree':   0.7,\n    'colsample_bylevel':  0.7,\n    'silent':             1\n}","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:41:41.889917Z","iopub.execute_input":"2022-03-29T14:41:41.890622Z","iopub.status.idle":"2022-03-29T14:41:41.897145Z","shell.execute_reply.started":"2022-03-29T14:41:41.890585Z","shell.execute_reply":"2022-03-29T14:41:41.896445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nxgbm = xgb.XGBRegressor(**xgb_params,num_boost_round = 200)\nxgbm.fit(X_train, y_train)\n# print(xgbm.score(X_train, y_train), xgbm.score(X_test, y_test), .score(test, trip_duration))\n","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:41:41.898444Z","iopub.execute_input":"2022-03-29T14:41:41.899229Z","iopub.status.idle":"2022-03-29T14:42:40.155298Z","shell.execute_reply.started":"2022-03-29T14:41:41.899192Z","shell.execute_reply":"2022-03-29T14:42:40.154596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ndtrain = xgb.DMatrix(X_train, label=y_train)\ndtest = xgb.DMatrix(X_test)\nparms = {'max_depth':15, #maximum depth of a tree 8 12\n         'objective':'reg:linear',\n         'eta'      :0.05, #0.3\n         'subsample':0.9,#SGD will use this percentage of data 0.8 0.99\n         'lambda'  :3, #L2 regularization term,>1 more conservative 4 \n         'colsample_bytree':0.6, #0.9\n         'colsample_bylevel':0.7, #1 0.7\n         'min_child_weight': 0.5, #10 0.5\n         #'nthread'  :3 ... default is max cores\n         'eval_metric':'rmse'}  #number of cpu core to use\n# running for 2k iterations \nxgbm2 = xgb.train(parms, dtrain, num_boost_round=125, maximize=False, verbose_eval=100)\nprint(np.sqrt(MSE(y_test, np.clip(xgbm2.predict(dtest),0,None))), end=\"\\n\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:42:40.156631Z","iopub.execute_input":"2022-03-29T14:42:40.157356Z","iopub.status.idle":"2022-03-29T14:43:56.511934Z","shell.execute_reply.started":"2022-03-29T14:42:40.157317Z","shell.execute_reply":"2022-03-29T14:43:56.511172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"models = [lgbm1,lgbm2,xgbm]","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:43:59.848222Z","iopub.execute_input":"2022-03-29T14:43:59.848489Z","iopub.status.idle":"2022-03-29T14:43:59.852107Z","shell.execute_reply.started":"2022-03-29T14:43:59.848451Z","shell.execute_reply":"2022-03-29T14:43:59.851439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for model in models:\n    print(str(model))\n    print(np.sqrt(MSE(y_train, np.clip(model.predict(X_train),0,None))), end=\"\\n\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:43:59.85331Z","iopub.execute_input":"2022-03-29T14:43:59.853746Z","iopub.status.idle":"2022-03-29T14:44:23.827301Z","shell.execute_reply.started":"2022-03-29T14:43:59.853708Z","shell.execute_reply":"2022-03-29T14:44:23.825808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# d_add_test = xgb.DMatrix(test)\nprint(np.sqrt(MSE(y_train, np.clip(xgbm2.predict(dtrain),0,None))), end=\"\\n\\n\")\n\nmodels.append(xgbm2)\nmodels","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:44:23.82871Z","iopub.execute_input":"2022-03-29T14:44:23.828962Z","iopub.status.idle":"2022-03-29T14:44:23.832406Z","shell.execute_reply.started":"2022-03-29T14:44:23.828926Z","shell.execute_reply":"2022-03-29T14:44:23.831704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Explanations","metadata":{}},{"cell_type":"code","source":"_, exp_data, _, y_exp_data = train_test_split(X, y, test_size=10000, random_state=42)\n","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:44:23.850454Z","iopub.execute_input":"2022-03-29T14:44:23.850691Z","iopub.status.idle":"2022-03-29T14:44:24.128656Z","shell.execute_reply.started":"2022-03-29T14:44:23.850657Z","shell.execute_reply":"2022-03-29T14:44:24.127809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Explainer","metadata":{}},{"cell_type":"code","source":"def get_explainer(model):\n    exp = shap.Explainer(model,exp_data)\n    return exp(exp_data,check_additivity=False)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:44:24.129989Z","iopub.execute_input":"2022-03-29T14:44:24.130228Z","iopub.status.idle":"2022-03-29T14:44:24.135086Z","shell.execute_reply.started":"2022-03-29T14:44:24.130195Z","shell.execute_reply":"2022-03-29T14:44:24.134391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lgbm1_exp = get_explainer(lgbm1)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:44:24.136358Z","iopub.execute_input":"2022-03-29T14:44:24.136765Z","iopub.status.idle":"2022-03-29T14:47:59.132626Z","shell.execute_reply.started":"2022-03-29T14:44:24.136724Z","shell.execute_reply":"2022-03-29T14:47:59.131901Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.plots.bar(lgbm1_exp)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:47:59.133909Z","iopub.execute_input":"2022-03-29T14:47:59.13423Z","iopub.status.idle":"2022-03-29T14:47:59.423146Z","shell.execute_reply.started":"2022-03-29T14:47:59.134192Z","shell.execute_reply":"2022-03-29T14:47:59.422495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lgbm2_exp = get_explainer(lgbm2)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T14:47:59.424362Z","iopub.execute_input":"2022-03-29T14:47:59.424779Z","iopub.status.idle":"2022-03-29T15:01:04.262735Z","shell.execute_reply.started":"2022-03-29T14:47:59.424741Z","shell.execute_reply":"2022-03-29T15:01:04.262014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.plots.bar(lgbm2_exp)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:01:04.264055Z","iopub.execute_input":"2022-03-29T15:01:04.264297Z","iopub.status.idle":"2022-03-29T15:01:04.539883Z","shell.execute_reply.started":"2022-03-29T15:01:04.264262Z","shell.execute_reply":"2022-03-29T15:01:04.539234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgbm_exp = get_explainer(xgbm)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:01:04.541194Z","iopub.execute_input":"2022-03-29T15:01:04.541677Z","iopub.status.idle":"2022-03-29T15:07:41.754027Z","shell.execute_reply.started":"2022-03-29T15:01:04.541631Z","shell.execute_reply":"2022-03-29T15:07:41.753209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.plots.bar(xgbm_exp)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:07:41.755312Z","iopub.execute_input":"2022-03-29T15:07:41.755577Z","iopub.status.idle":"2022-03-29T15:07:42.026144Z","shell.execute_reply.started":"2022-03-29T15:07:41.755542Z","shell.execute_reply":"2022-03-29T15:07:42.025489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgbm2_exp = get_explainer(xgbm2)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:07:42.027177Z","iopub.execute_input":"2022-03-29T15:07:42.027577Z","iopub.status.idle":"2022-03-29T15:17:39.560364Z","shell.execute_reply.started":"2022-03-29T15:07:42.027535Z","shell.execute_reply":"2022-03-29T15:17:39.559636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.plots.bar(xgbm2_exp)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:17:39.561759Z","iopub.execute_input":"2022-03-29T15:17:39.562005Z","iopub.status.idle":"2022-03-29T15:17:39.833954Z","shell.execute_reply.started":"2022-03-29T15:17:39.56197Z","shell.execute_reply":"2022-03-29T15:17:39.833277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Shap values","metadata":{}},{"cell_type":"code","source":"def get_shap_vals(model):\n    exp = shap.Explainer(model,X_train)\n    shap_vals = exp.shap_values(exp_data,check_additivity=False)\n    return shap_vals","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:17:39.835058Z","iopub.execute_input":"2022-03-29T15:17:39.835372Z","iopub.status.idle":"2022-03-29T15:17:39.839673Z","shell.execute_reply.started":"2022-03-29T15:17:39.835335Z","shell.execute_reply":"2022-03-29T15:17:39.838957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap_vals = np.empty((len(models), *exp_data.shape))\nfor i,model in enumerate(models):\n    print(i)\n    shap_vals[i] = get_shap_vals(model)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:17:39.85104Z","iopub.execute_input":"2022-03-29T15:17:39.851716Z","iopub.status.idle":"2022-03-29T15:50:24.731598Z","shell.execute_reply.started":"2022-03-29T15:17:39.851677Z","shell.execute_reply":"2022-03-29T15:50:24.730704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap_vals.shape","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:24.732888Z","iopub.execute_input":"2022-03-29T15:50:24.73314Z","iopub.status.idle":"2022-03-29T15:50:24.739092Z","shell.execute_reply.started":"2022-03-29T15:50:24.733107Z","shell.execute_reply":"2022-03-29T15:50:24.738454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_shap_vals = np.empty((4, 10000,16))","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:24.757825Z","iopub.execute_input":"2022-03-29T15:50:24.758063Z","iopub.status.idle":"2022-03-29T15:50:24.767421Z","shell.execute_reply.started":"2022-03-29T15:50:24.75803Z","shell.execute_reply":"2022-03-29T15:50:24.76673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_shap_vals[0] = np.genfromtxt(\"../input/nyctaxishapvals/lgbm1_shap_vals.csv\", delimiter=',')  \nall_shap_vals[1] = np.genfromtxt(\"../input/nyctaxishapvals/lgbm2_shap_vals.csv\", delimiter=',')  \nall_shap_vals[2] = np.genfromtxt(\"../input/nyctaxishapvals/xgbm1_shap_vals.csv\", delimiter=',')  \nall_shap_vals[3] = np.genfromtxt(\"../input/nyctaxishapvals/xgbm2_shap_vals.csv\", delimiter=',')  ","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:24.770019Z","iopub.execute_input":"2022-03-29T15:50:24.770254Z","iopub.status.idle":"2022-03-29T15:50:25.577991Z","shell.execute_reply.started":"2022-03-29T15:50:24.770224Z","shell.execute_reply":"2022-03-29T15:50:25.577271Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_shap_vals.shape","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:25.579131Z","iopub.execute_input":"2022-03-29T15:50:25.579403Z","iopub.status.idle":"2022-03-29T15:50:25.58537Z","shell.execute_reply.started":"2022-03-29T15:50:25.57936Z","shell.execute_reply":"2022-03-29T15:50:25.584561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### similarity","metadata":{}},{"cell_type":"code","source":"def shap_similarity(epsilon=20, normalize=True, measure=0):\n    sim_matrix = np.zeros([len(models),len(models)])\n    for i in range(len(models)-1):\n        for j in range(i+1,len(models)):\n            shap_vals1 = shap_vals[i]\n            shap_vals2 = shap_vals[j]\n            if measure == 0:\n                similarity = np.count_nonzero(np.abs(shap_vals1-shap_vals2) < epsilon) \n                if normalize:\n                    similarity /= np.prod(shap_vals1.shape) \n            elif measure == 1:\n                similarity = np.count_nonzero(np.sum(np.abs(shap_vals1-shap_vals2),axis=1)<epsilon*shap_vals1.shape[1])\n                if normalize:\n                    similarity /= shap_vals1.shape[0]\n            elif measure == \"loss\":\n                similarity = np.sum(np.abs(shap_vals1-shap_vals2))\n                if normalize:\n                    similarity /= np.prod(shap_vals1.shape) \n            else:\n                similarity = np.count_nonzero(np.max(np.abs(shap_vals1-shap_vals2),axis=1)<epsilon)\n                if normalize:\n                    similarity /= shap_vals1.shape[0]\n            print(str(models[i])[:120],\"vs\",str(models[j])[:120],similarity, sep=\"\\n\", end=\"\\n\\n\")\n            sim_matrix[i,j] = sim_matrix[j, i] = similarity    \n    return np.sum(sim_matrix,axis=0) / (len(models)-1)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:25.586685Z","iopub.execute_input":"2022-03-29T15:50:25.58717Z","iopub.status.idle":"2022-03-29T15:50:25.600246Z","shell.execute_reply.started":"2022-03-29T15:50:25.587131Z","shell.execute_reply":"2022-03-29T15:50:25.599395Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap_similarity()","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:25.601319Z","iopub.execute_input":"2022-03-29T15:50:25.601752Z","iopub.status.idle":"2022-03-29T15:50:25.62785Z","shell.execute_reply.started":"2022-03-29T15:50:25.601716Z","shell.execute_reply":"2022-03-29T15:50:25.627129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"delta = [10, 12.5, 15, 17.5, 20]","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:25.62891Z","iopub.execute_input":"2022-03-29T15:50:25.629519Z","iopub.status.idle":"2022-03-29T15:50:25.633643Z","shell.execute_reply.started":"2022-03-29T15:50:25.629484Z","shell.execute_reply":"2022-03-29T15:50:25.632702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#T-sim \nT_delta_20 =  [20, 0.7546    , 0.69835   , 0.76955833, 0.7540625]\nT_delta_175 = [17.5, 0.72081458, 0.6654    , 0.73611042, 0.7203  ]\nT_delta_15 =  [15, 0.68227292, 0.62793125, 0.69672292, 0.68137292]\nT_delta_125 = [12.5, 0.6366875 , 0.58493958, 0.65024583, 0.63460625]\nT_delta_10 =  [10, 0.58339583, 0.53530833, 0.59498125, 0.58037708]\nrows = [T_delta_20 ,\n        T_delta_175,\n        T_delta_15 ,\n        T_delta_125,\n        T_delta_10 ,]\ncols = [\"delta\", \"lgbm1\", \"lgbm2\" ,\"xgbm1\",\"xgbm2\"]\nT_delta_df = pd.DataFrame(rows,columns=cols)\nT_delta_df","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:25.634973Z","iopub.execute_input":"2022-03-29T15:50:25.635506Z","iopub.status.idle":"2022-03-29T15:50:25.651539Z","shell.execute_reply.started":"2022-03-29T15:50:25.635472Z","shell.execute_reply":"2022-03-29T15:50:25.650794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\nsns.set(rc = {'figure.figsize':(15,8)})\nsns.set_theme(style=\"whitegrid\")\n    ","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:25.652673Z","iopub.execute_input":"2022-03-29T15:50:25.652981Z","iopub.status.idle":"2022-03-29T15:50:25.706025Z","shell.execute_reply.started":"2022-03-29T15:50:25.652947Z","shell.execute_reply":"2022-03-29T15:50:25.705366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sim_fig = sns.lineplot(x=\"delta\", y=\"value\", hue=\"variable\", style =\"variable\", markers=True, data=pd.melt(T_delta_df,[\"delta\"]))\nsim_fig.set_ylabel(\"average similarity coefficient\")\nsim_fig.set_yticks(np.linspace(0.5,0.8,11))","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:25.707062Z","iopub.execute_input":"2022-03-29T15:50:25.707383Z","iopub.status.idle":"2022-03-29T15:50:26.08776Z","shell.execute_reply.started":"2022-03-29T15:50:25.707339Z","shell.execute_reply":"2022-03-29T15:50:26.087032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#S-simlarity \ndelta_10 =  [10,0.24066667, 0.07083333, 0.3572    , 0.3327   ]\ndelta_125 = [12.5, 0.46983333, 0.21253333, 0.53946667, 0.4827 ]\ndelta_15 =  [15, 0.64136667, 0.40136667, 0.67413333, 0.59986667]\ndelta_175 = [17.5, 0.75196667, 0.579     , 0.77223333, 0.70686667]\ndelta_20 =  [20, 0.81503333, 0.69933333, 0.83596667, 0.792 ]\n\nrows = [delta_20 ,\n        delta_175,\n        delta_15 ,\n        delta_125,\n        delta_10 ,]\ncols = [\"delta\", \"lgbm1\", \"lgbm2\" ,\"xgbm1\",\"xgbm2\"]\ndelta_df = pd.DataFrame(rows,columns=cols)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.088747Z","iopub.execute_input":"2022-03-29T15:50:26.088993Z","iopub.status.idle":"2022-03-29T15:50:26.096498Z","shell.execute_reply.started":"2022-03-29T15:50:26.088956Z","shell.execute_reply":"2022-03-29T15:50:26.095638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sim_fig = sns.lineplot(x=\"delta\", y=\"value\", hue=\"variable\", style =\"variable\", markers=True, data=pd.melt(delta_df,[\"delta\"]))\nsim_fig.set_ylabel(\"average similarity coefficient\")","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.097833Z","iopub.execute_input":"2022-03-29T15:50:26.098281Z","iopub.status.idle":"2022-03-29T15:50:26.485756Z","shell.execute_reply.started":"2022-03-29T15:50:26.098235Z","shell.execute_reply":"2022-03-29T15:50:26.485077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"####  conflict","metadata":{}},{"cell_type":"code","source":"def shap_conflict(alpha=400,normalize=True,measure=\"gran\"):\n    conf_matrix = np.zeros([len(models),len(models)])\n    for i in range(len(models)-1):\n        for j in range(i+1,len(models)):\n            shap_vals1 = shap_vals[i]\n            shap_vals2 = shap_vals[j]\n            conf_cond = shap_vals1*shap_vals2<0\n            dist_cond = np.abs(shap_vals1-shap_vals2)> alpha\n            satifies_conditions = np.logical_and(dist_cond,conf_cond)\n            print(np.nonzero(satifies_conditions))\n            if measure == \"gran\":\n                conflict = np.sum(satifies_conditions)\n                if normalize:\n                    conflict /= np.prod(shap_vals1.shape)\n            elif measure == \"def\":\n                conflict = np.count_nonzero(np.max(satifies_conditions, axis=1))\n                if normalize:\n                    conflict /= shap_vals1.shape[0]\n            print(str(models[i])[:120],\"vs\",str(models[j])[:120],conflict, sep=\"\\n\", end=\"\\n\\n\")\n            conf_matrix[i,j] = conf_matrix[j, i] = conflict    \n    return np.sum(conf_matrix,axis=0) / (len(models)-1)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.48713Z","iopub.execute_input":"2022-03-29T15:50:26.487383Z","iopub.status.idle":"2022-03-29T15:50:26.496788Z","shell.execute_reply.started":"2022-03-29T15:50:26.487345Z","shell.execute_reply":"2022-03-29T15:50:26.496014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap_conflict()","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.498307Z","iopub.execute_input":"2022-03-29T15:50:26.498584Z","iopub.status.idle":"2022-03-29T15:50:26.534775Z","shell.execute_reply.started":"2022-03-29T15:50:26.498548Z","shell.execute_reply":"2022-03-29T15:50:26.534129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap_vals[0][5725,3]","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.536277Z","iopub.execute_input":"2022-03-29T15:50:26.536528Z","iopub.status.idle":"2022-03-29T15:50:26.541877Z","shell.execute_reply.started":"2022-03-29T15:50:26.536493Z","shell.execute_reply":"2022-03-29T15:50:26.541129Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap_vals[1][5725,3]","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.543364Z","iopub.execute_input":"2022-03-29T15:50:26.54395Z","iopub.status.idle":"2022-03-29T15:50:26.551303Z","shell.execute_reply.started":"2022-03-29T15:50:26.543915Z","shell.execute_reply":"2022-03-29T15:50:26.550442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_index = 1184","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.552833Z","iopub.execute_input":"2022-03-29T15:50:26.553136Z","iopub.status.idle":"2022-03-29T15:50:26.557576Z","shell.execute_reply.started":"2022-03-29T15:50:26.553101Z","shell.execute_reply":"2022-03-29T15:50:26.556762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.bar_plot(shap_vals[0][data_index],X.iloc[[data_index]])","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.565702Z","iopub.execute_input":"2022-03-29T15:50:26.565889Z","iopub.status.idle":"2022-03-29T15:50:26.772927Z","shell.execute_reply.started":"2022-03-29T15:50:26.565866Z","shell.execute_reply":"2022-03-29T15:50:26.772278Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.bar_plot(shap_vals[1][data_index],X.iloc[[data_index]])","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.774095Z","iopub.execute_input":"2022-03-29T15:50:26.774787Z","iopub.status.idle":"2022-03-29T15:50:26.967972Z","shell.execute_reply.started":"2022-03-29T15:50:26.774748Z","shell.execute_reply":"2022-03-29T15:50:26.967313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shap.summary_plot(shap_vals[0], exp_data)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:26.969239Z","iopub.execute_input":"2022-03-29T15:50:26.969509Z","iopub.status.idle":"2022-03-29T15:50:29.045868Z","shell.execute_reply.started":"2022-03-29T15:50:26.969473Z","shell.execute_reply":"2022-03-29T15:50:29.04519Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"alpha = [80,100,120,140,160]","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:29.047081Z","iopub.execute_input":"2022-03-29T15:50:29.047452Z","iopub.status.idle":"2022-03-29T15:50:29.052395Z","shell.execute_reply.started":"2022-03-29T15:50:29.047399Z","shell.execute_reply":"2022-03-29T15:50:29.051671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#C-conf\nalpha_80 = np.array([80,0.06383333, 0.1183    , 0.0439    , 0.05123333])\nalpha_100 =np.array([100,0.04486667, 0.0859    , 0.02773333, 0.0335    ])\nalpha_120= np.array([120,0.03376667, 0.06806667, 0.01966667, 0.02516667])\nalpha_140 =np.array([140,0.0262    , 0.05486667, 0.0146    , 0.0196    ])\nalpha_160 =np.array([160,0.0219    , 0.04466667, 0.0111    , 0.01473333])\n\n\nrows = [alpha_80 ,\n        alpha_100,\n        alpha_120,\n        alpha_140,\n        alpha_160,]\ncols = [\"alpha\", \"lgbm1\", \"lgbm2\" ,\"xgbm1\",\"xgbm2\"]\nalpha_df = pd.DataFrame(rows,columns=cols)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:29.053734Z","iopub.execute_input":"2022-03-29T15:50:29.054544Z","iopub.status.idle":"2022-03-29T15:50:29.063694Z","shell.execute_reply.started":"2022-03-29T15:50:29.05451Z","shell.execute_reply":"2022-03-29T15:50:29.062941Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"C_conf_fig = sns.lineplot(x=\"alpha\", y=\"value\", hue=\"variable\", style =\"variable\", markers=True, data=pd.melt(alpha_df,[\"alpha\"]))\nC_conf_fig.set_ylabel(\"average conflict coefficient\")\nC_conf_fig.set_yticks(np.linspace(0.01,0.12,12))","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:29.064834Z","iopub.execute_input":"2022-03-29T15:50:29.065162Z","iopub.status.idle":"2022-03-29T15:50:29.641136Z","shell.execute_reply.started":"2022-03-29T15:50:29.065126Z","shell.execute_reply":"2022-03-29T15:50:29.640381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#F-conf\nF_alpha_80  = np.array([80,0.00433333, 0.00865208, 0.00318542, 0.00385417])\nF_alpha_100 = np.array([100,0.00297292, 0.00611875, 0.00195625, 0.00248958])\nF_alpha_120 = np.array([120,0.00222708, 0.0047125 , 0.00135833, 0.00180625])\nF_alpha_140 = np.array([140,0.00172292, 0.00374167, 0.00098333, 0.00138542])\nF_alpha_160 = np.array([160,0.00143542, 0.00300833, 0.00074167, 0.00102292])\n\nrows = [F_alpha_80 ,\n        F_alpha_100,\n        F_alpha_120,\n        F_alpha_140,\n        F_alpha_160,]\ncols = [\"alpha\", \"lgbm1\", \"lgbm2\" ,\"xgbm1\",\"xgbm2\"]\nF_alpha_df = pd.DataFrame(rows,columns=cols)","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:29.642121Z","iopub.execute_input":"2022-03-29T15:50:29.642369Z","iopub.status.idle":"2022-03-29T15:50:29.656891Z","shell.execute_reply.started":"2022-03-29T15:50:29.642334Z","shell.execute_reply":"2022-03-29T15:50:29.655957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"F_conf_fig = sns.lineplot(x=\"alpha\", y=\"value\", hue=\"variable\", style =\"variable\", markers=True, data=pd.melt(F_alpha_df,[\"alpha\"]))\nF_conf_fig.set_ylabel(\"average conflict coefficient\")","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:29.659096Z","iopub.execute_input":"2022-03-29T15:50:29.659333Z","iopub.status.idle":"2022-03-29T15:50:30.089965Z","shell.execute_reply.started":"2022-03-29T15:50:29.6593Z","shell.execute_reply":"2022-03-29T15:50:30.089289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generalisation  hypothesis","metadata":{}},{"cell_type":"code","source":"best_sim_model = #\nworst_sim_model = #","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:30.091161Z","iopub.execute_input":"2022-03-29T15:50:30.093016Z","iopub.status.idle":"2022-03-29T15:50:30.100394Z","shell.execute_reply.started":"2022-03-29T15:50:30.092987Z","shell.execute_reply":"2022-03-29T15:50:30.098717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(str(best_sim_model))\nprint(np.sqrt(MSE(y_extra, np.clip(best_sim_model.predict(X_extra),0,None))), end=\"\\n\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:30.101261Z","iopub.status.idle":"2022-03-29T15:50:30.102081Z","shell.execute_reply.started":"2022-03-29T15:50:30.10188Z","shell.execute_reply":"2022-03-29T15:50:30.101907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(str(worst_sim_model))\nprint(np.sqrt(MSE(y_extra, np.clip(worst_sim_model.predict(X_extra),0,None))), end=\"\\n\\n\")","metadata":{"execution":{"iopub.status.busy":"2022-03-29T15:50:30.103624Z","iopub.status.idle":"2022-03-29T15:50:30.104259Z","shell.execute_reply.started":"2022-03-29T15:50:30.104024Z","shell.execute_reply":"2022-03-29T15:50:30.104051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### d. Cross-validation <a id=\"four-d\"></a>","metadata":{"_uuid":"122dcd986e18d754077660fff653c542231d6a36"}},{"cell_type":"code","source":"# lgb_df = lgb.Dataset(X, y)\n# lgb.cv(lgb_params, lgb_df, stratified=False) #False is needed as it only works with classification\n\n# # Cross-validation on LightGBM model (sklearn API) ------------\n# from sklearn.model_selection import cross_val_score\n\n# lgb1_cv_score = cross_val_score(lgbm, X, y, cv=5, scoring=\"neg_mean_squared_log_error\")\n# print(lgb_cv_score)\n# print(np.mean(lgb_cv_score))\n\n# xgb_cv_score = cross_val_score(xgbm, X, y, cv=5, scoring=\"neg_mean_squared_log_error\")\n# print(xgb_cv_score)\n# print(np.mean(xgb_cv_score))\n\n# #Output:\n#     #[0.77872018 0.7801329  0.77988107 0.78049745 0.77904688]\n#     #0.7796556968369478","metadata":{"_uuid":"57ec31cb4ee2cf52963732cb9337818478f6c3b0","execution":{"iopub.status.busy":"2022-03-29T15:50:30.105417Z","iopub.status.idle":"2022-03-29T15:50:30.105999Z","shell.execute_reply.started":"2022-03-29T15:50:30.105761Z","shell.execute_reply":"2022-03-29T15:50:30.105788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# SRA results","metadata":{}},{"cell_type":"code","source":"def top_features(datapoint_shap_vals, nfeatures=7):\n    \"\"\"gives the indexes of the n features with the biggest absolute shap vals\"\"\"\n    return np.argpartition(abs(datapoint_shap_vals), -nfeatures)[-nfeatures:]\n    \n\ndef equal_ranking(feat, top_feats1, top_feats2):\n    return list(np.where(top_feats1==feat)[0]) == list(np.where(top_feats2==feat)[0])\n\ndef aggregated_sra(k=3, agg=np.mean):\n    sra_matrix = np.zeros([len(all_shap_vals),len(all_shap_vals)])\n    for i in range(len(all_shap_vals)-1):\n        for j in range(i+1,len(all_shap_vals)):\n            shap_vals1 = all_shap_vals[i]\n            shap_vals2 = all_shap_vals[j]\n            ndatapoints = len(shap_vals1)\n            model_sras = np.empty([ndatapoints])\n            for l in range(ndatapoints):\n                sra = 0\n                top_feats1 = top_features(shap_vals1[l], nfeatures=k)\n                top_feats2 = top_features(shap_vals2[l], nfeatures=k)\n                for feat in top_feats1:\n                    if equal_ranking(feat, top_feats1, top_feats2) and shap_vals1[l][feat]*shap_vals2[l][feat]>=0:\n                        sra += 1\n                sra /= k\n                model_sras[l] = sra\n            sra_matrix[i, j] = sra_matrix[j,i] =agg(model_sras)\n    print(sra_matrix)\n    return np.sum(sra_matrix,axis=0) / (len(all_shap_vals)-1)\n","metadata":{"execution":{"iopub.status.busy":"2022-08-11T15:50:29.072104Z","iopub.execute_input":"2022-08-11T15:50:29.073098Z","iopub.status.idle":"2022-08-11T15:50:29.163995Z","shell.execute_reply.started":"2022-08-11T15:50:29.072998Z","shell.execute_reply":"2022-08-11T15:50:29.162676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"aggregated_sra()","metadata":{},"execution_count":null,"outputs":[]}]}