{"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":"# Imports","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport pyproj\ngeodesic = pyproj.Geod(ellps='WGS84')\n\nfrom sklearn.preprocessing import OrdinalEncoder, PolynomialFeatures\nfrom sklearn.compose import ColumnTransformer\nfrom sklearn.model_selection import train_test_split, GridSearchCV\nimport xgboost as xgb","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:05.147752Z","iopub.execute_input":"2022-08-10T17:00:05.148116Z","iopub.status.idle":"2022-08-10T17:00:05.156657Z","shell.execute_reply.started":"2022-08-10T17:00:05.148085Z","shell.execute_reply":"2022-08-10T17:00:05.155126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nrows = 5*10**6\ntrain = pd.read_csv('../input/new-york-city-taxi-fare-prediction/train.csv', nrows = nrows)\ntest  = pd.read_csv('../input/new-york-city-taxi-fare-prediction/test.csv')","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:05.158536Z","iopub.execute_input":"2022-08-10T17:00:05.159044Z","iopub.status.idle":"2022-08-10T17:00:17.685142Z","shell.execute_reply.started":"2022-08-10T17:00:05.158999Z","shell.execute_reply":"2022-08-10T17:00:17.683531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Glimpse the data and deal with NA values","metadata":{}},{"cell_type":"code","source":"train.info(show_counts = True)\nprint('-'*50)\ntest.info()","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:17.686921Z","iopub.execute_input":"2022-08-10T17:00:17.687433Z","iopub.status.idle":"2022-08-10T17:00:18.211052Z","shell.execute_reply.started":"2022-08-10T17:00:17.687383Z","shell.execute_reply":"2022-08-10T17:00:18.210057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* There are no NA values in the test set. So, NA values from the train set can be considered as unimportant and can be dropped.","metadata":{}},{"cell_type":"code","source":"train = train.dropna()","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:18.212386Z","iopub.execute_input":"2022-08-10T17:00:18.213631Z","iopub.status.idle":"2022-08-10T17:00:19.279469Z","shell.execute_reply.started":"2022-08-10T17:00:18.213577Z","shell.execute_reply":"2022-08-10T17:00:19.278018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.describe().round(3)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:19.281026Z","iopub.execute_input":"2022-08-10T17:00:19.281404Z","iopub.status.idle":"2022-08-10T17:00:20.722750Z","shell.execute_reply.started":"2022-08-10T17:00:19.281371Z","shell.execute_reply":"2022-08-10T17:00:20.721472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test.describe().round(3)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:20.724391Z","iopub.execute_input":"2022-08-10T17:00:20.724790Z","iopub.status.idle":"2022-08-10T17:00:20.763171Z","shell.execute_reply.started":"2022-08-10T17:00:20.724757Z","shell.execute_reply":"2022-08-10T17:00:20.761648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* minimum value of fare is negative. It might be data recording error. I will check it in more depth after calculating the distance between coordinates.\n* In the training set there are outliers in coordinate features, but the test data doesn't suffer from this.\n* Description shows that minimum and maximum number of passengers have unrealistic values in the training set.","metadata":{}},{"cell_type":"markdown","source":"# Passenger count and coordinate features","metadata":{}},{"cell_type":"markdown","source":"Let's make range of coordinate features from train and test set equal to each other($\\pm$ 0.5 lat/long). and then explore passenger_count further.","metadata":{}},{"cell_type":"code","source":"coordinate_features = ['pickup_longitude', 'pickup_latitude', 'dropoff_longitude', 'dropoff_latitude']\n\nfor feature in coordinate_features:\n    minn = test[feature].min()\n    maxx = test[feature].max()\n    stdd = test[feature].std()\n    \n    train = train[((minn - 0.5) < train[feature]) & (train[feature] < (maxx + 0.5))]","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:20.764910Z","iopub.execute_input":"2022-08-10T17:00:20.765356Z","iopub.status.idle":"2022-08-10T17:00:24.318303Z","shell.execute_reply.started":"2022-08-10T17:00:20.765320Z","shell.execute_reply":"2022-08-10T17:00:24.314815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.style.use('seaborn-whitegrid')\nfig, axes = plt.subplots(nrows=1,ncols=2, figsize = (18,5))\ntrain['passenger_count'].value_counts(normalize = 'columns').plot(kind = 'bar', ax = axes[0])\ntest['passenger_count'].value_counts(normalize = 'columns').plot(kind = 'bar', ax = axes[1])\naxes[0].set_title('Distribution of passengers [Train]', fontsize = 15)\naxes[1].set_title('Distribution of passengers [Test]', fontsize = 15);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:24.323172Z","iopub.execute_input":"2022-08-10T17:00:24.325791Z","iopub.status.idle":"2022-08-10T17:00:24.992315Z","shell.execute_reply.started":"2022-08-10T17:00:24.325653Z","shell.execute_reply":"2022-08-10T17:00:24.990883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Seems data recording errors. All the values which are only present in the train data but not in the test data can be removed safely.","metadata":{}},{"cell_type":"code","source":"train = train[train['passenger_count'].isin(test['passenger_count'].unique())] # Removing\n\npd.pivot_table(train, 'fare_amount', 'passenger_count', aggfunc = ['mean', 'std'])","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:24.996091Z","iopub.execute_input":"2022-08-10T17:00:24.996800Z","iopub.status.idle":"2022-08-10T17:00:26.305383Z","shell.execute_reply.started":"2022-08-10T17:00:24.996762Z","shell.execute_reply":"2022-08-10T17:00:26.303623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* passenger_count seems redundant for our prediction.","metadata":{}},{"cell_type":"code","source":"print(f'So far we have lost {round(1-train.shape[0]/nrows,4)*100}% of observations from the training set.')","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:26.307367Z","iopub.execute_input":"2022-08-10T17:00:26.307780Z","iopub.status.idle":"2022-08-10T17:00:26.316199Z","shell.execute_reply.started":"2022-08-10T17:00:26.307746Z","shell.execute_reply":"2022-08-10T17:00:26.314283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Distance and azimuth","metadata":{}},{"cell_type":"markdown","source":"The obvious thing to do from two coordinates is to calculate the distance and the azimuth between them.","metadata":{}},{"cell_type":"code","source":"test['azimuth'], cache, test['travelled_m'] = geodesic.inv(\n    test['pickup_latitude'].values, test['pickup_longitude'].values,\n    test['dropoff_latitude'].values, test['dropoff_longitude'].values\n)\ntrain['azimuth'], cache, train['travelled_m'] = geodesic.inv(\n    train['pickup_latitude'].values, train['pickup_longitude'].values,\n    train['dropoff_latitude'].values, train['dropoff_longitude'].values\n)\n\n#From meters to kilometers\ntrain['travelled_km'] = train['travelled_m']/1000 \ntest['travelled_km'] = test['travelled_m']/1000 ","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:00:58.921618Z","iopub.execute_input":"2022-08-10T17:00:58.923817Z","iopub.status.idle":"2022-08-10T17:01:04.466181Z","shell.execute_reply.started":"2022-08-10T17:00:58.923728Z","shell.execute_reply":"2022-08-10T17:01:04.464576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(figsize=(15,4), nrows = 1, ncols = 2)\naxes[0].set_title('Distribution of Travelled km [Train]', fontsize = 13)\naxes[0].set_xlabel('travelled_km')\naxes[0].set_ylabel('Density')\naxes[0].hist(train['travelled_km'], bins = 1000, density = True);\n\naxes[1].set_title('Distribution of Travelled km [Test]', fontsize = 13)\naxes[1].set_xlabel('travelled_km')\naxes[1].set_ylabel('Density')\naxes[1].hist(test['travelled_km'], bins = 1000, density = True);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:01:07.292527Z","iopub.execute_input":"2022-08-10T17:01:07.293301Z","iopub.status.idle":"2022-08-10T17:01:11.710110Z","shell.execute_reply.started":"2022-08-10T17:01:07.293248Z","shell.execute_reply":"2022-08-10T17:01:11.708609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* There are definitely outliers in both datasets.\n* Exploring the relation between distance and fare will give us a better idea of how to deal with outliers.","metadata":{}},{"cell_type":"code","source":"plt.scatter(train['travelled_km'], train['fare_amount'], alpha = 0.1)\nplt.xlabel('travelled_km')\nplt.ylabel('Fare amount')\nplt.title('Scatterplot of Fare amount vs travelled_km');","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:02:05.007748Z","iopub.execute_input":"2022-08-10T17:02:05.008188Z","iopub.status.idle":"2022-08-10T17:02:14.369980Z","shell.execute_reply.started":"2022-08-10T17:02:05.008155Z","shell.execute_reply":"2022-08-10T17:02:14.368496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Observations, where distances are more than 70km, will be replaced by the median distance, as most of them seem to be data recording errors.\n* For stability of our model, observations, where fares are more than 300, will be removed.\n* From the graph, it's visible that some people paid a negative or zero amount which is unrealistic. Let's zoom into those values.","metadata":{}},{"cell_type":"code","source":"train.loc[train['travelled_km'] > 70, 'travelled_km'] = train['travelled_km'].median() \ntest.loc[test['travelled_km'] > 70, 'travelled_km']   = test['travelled_km'].median()\ntrain = train[train['fare_amount'] < 300]","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:02:26.177707Z","iopub.execute_input":"2022-08-10T17:02:26.178185Z","iopub.status.idle":"2022-08-10T17:02:27.323030Z","shell.execute_reply.started":"2022-08-10T17:02:26.178146Z","shell.execute_reply":"2022-08-10T17:02:27.321838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.scatter(\n    train.loc[train['fare_amount'] <= 0, 'travelled_km'],\n    train.loc[train['fare_amount'] <= 0, 'fare_amount']\n)\nplt.xlabel('Distance (km)')\nplt.ylabel('Fare amount')\nplt.title('Scatterplot for neagtive fare amount');","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2022-08-10T17:03:10.128724Z","iopub.execute_input":"2022-08-10T17:03:10.129944Z","iopub.status.idle":"2022-08-10T17:03:10.352309Z","shell.execute_reply.started":"2022-08-10T17:03:10.129891Z","shell.execute_reply":"2022-08-10T17:03:10.350889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* negative sign of fares seems to be just data recording errors. They will be changed by plus signs.\n* Fares equal to zeros or even less than ones seem unrealistic. They will be dropped.\n* There are some people who travelled 0 km but paid a positive amount. We have the same problem in the test set. So, I will leave 0 distances unchanged.","metadata":{}},{"cell_type":"code","source":"train['fare_amount'] = abs(train['fare_amount'])\ntrain = train[train['fare_amount'] >=1]","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:04:14.244276Z","iopub.execute_input":"2022-08-10T17:04:14.244812Z","iopub.status.idle":"2022-08-10T17:04:14.948379Z","shell.execute_reply.started":"2022-08-10T17:04:14.244767Z","shell.execute_reply":"2022-08-10T17:04:14.946478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axes = plt.subplots(figsize=(15,4), nrows=1, ncols=2)\naxes[0].set_title('Distribution of azimuth [Train]', fontsize = 13)\naxes[0].set_xlabel('azimuth')\naxes[0].set_ylabel('Density')\naxes[0].hist(train['azimuth'], bins = 100, density = True);\n\naxes[1].set_title('Distribution of azimuth [Test]', fontsize = 13)\naxes[1].set_xlabel('azimuth')\naxes[1].set_ylabel('Density')\naxes[1].hist(test['azimuth'], bins = 100, density = True);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:04:20.177443Z","iopub.execute_input":"2022-08-10T17:04:20.177871Z","iopub.status.idle":"2022-08-10T17:04:20.934350Z","shell.execute_reply.started":"2022-08-10T17:04:20.177839Z","shell.execute_reply":"2022-08-10T17:04:20.933035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Passengers are travelling mostly towards north or south (Maybe due to some geographical peculiarities of New York)","metadata":{}},{"cell_type":"code","source":"print(f'So far we have lost {round(1-train.shape[0]/nrows,4)*100}% of observations from the training set.')","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:04:25.304064Z","iopub.execute_input":"2022-08-10T17:04:25.304501Z","iopub.status.idle":"2022-08-10T17:04:25.311589Z","shell.execute_reply.started":"2022-08-10T17:04:25.304466Z","shell.execute_reply":"2022-08-10T17:04:25.310164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploring patterns on the map and adding new features","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(nrows=2,ncols=2, figsize = (24,20))\nmap_box = [-74.743, -72.509, 40.079, 42.208] # taken from the map.\nnyc_img = plt.imread('../input/nyc-map/map.png')\n\ntrain.plot(kind = 'scatter', x = 'pickup_longitude', y = 'pickup_latitude',\n           c = 'fare_amount', cmap = plt.get_cmap('jet'), colorbar = True,\n           alpha = 0.5, ax = axes[0,0])\n\naxes[0,0].set_title('Pickup locations', fontsize = 23)\naxes[0,0].set_xlabel('')\naxes[0,0].set_ylabel('')\naxes[0,0].imshow(nyc_img, extent = map_box, alpha = 1, cmap = plt.get_cmap('jet'), vmin = 0, vmax = 300)\n\ntrain[train['fare_amount'] > 80].plot(kind = 'scatter', x = 'pickup_longitude', y = 'pickup_latitude',\n                                      c = 'fare_amount', cmap = plt.get_cmap('jet'), colorbar = True,\n                                      alpha = 0.5, ax = axes[0,1])\n\naxes[0,1].set_title('Pickup locations[fare_amount > 80]', fontsize = 23)\naxes[0,1].set_xlabel('')\naxes[0,1].set_ylabel('')\naxes[0,1].imshow(nyc_img, extent = map_box, alpha=1, cmap=plt.get_cmap(\"jet\"), vmin = 0, vmax = 300)\n\ntrain.plot(kind = 'scatter', x = 'dropoff_longitude', y = 'dropoff_latitude',\n           c = 'fare_amount',cmap = plt.get_cmap('jet'), colorbar = True,\n           alpha = 0.5, ax = axes[1,0])\n\naxes[1,0].set_title('Dropoff locations', fontsize = 23)\naxes[1,0].set_xlabel('')\naxes[1,0].set_ylabel('')\naxes[1,0].imshow(nyc_img, extent = map_box, alpha = 1, cmap = plt.get_cmap('jet'), vmin = 0, vmax = 300)\n\ntrain[train['fare_amount'] > 80].plot(kind = 'scatter', x = 'dropoff_longitude', y = 'dropoff_latitude',\n                                      c = 'fare_amount', cmap = plt.get_cmap('jet'), colorbar = True,\n                                      alpha = 0.5, ax = axes[1,1])\n\naxes[1,1].set_title('Dropoff locations[fare_amount > 80]', fontsize = 23)\naxes[1,1].set_xlabel('')\naxes[1,1].set_ylabel('')\naxes[1,1].imshow(nyc_img, extent = map_box, alpha=1, cmap=plt.get_cmap('jet'), vmin = 0, vmax = 300);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T18:20:34.394001Z","iopub.execute_input":"2022-08-10T18:20:34.394434Z","iopub.status.idle":"2022-08-10T18:23:51.269464Z","shell.execute_reply.started":"2022-08-10T18:20:34.394399Z","shell.execute_reply":"2022-08-10T18:23:51.267769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The above image and coordinates are taken from [openstreetmap.org](https://www.openstreetmap.org/).\n* It seems there are some patterns on the map. Travellers picked up or dropped off far from the City of New York seem to have higher fares. \n\nLet's explore patterns further.","metadata":{}},{"cell_type":"code","source":"n_cat = 30\ntrain['pickup_longitude_cat']  = pd.cut(train['pickup_longitude'],  n_cat)\ntrain['pickup_latitude_cat']   = pd.cut(train['pickup_latitude'],   n_cat)\ntrain['dropoff_longitude_cat'] = pd.cut(train['dropoff_longitude'], n_cat)\ntrain['dropoff_latitude_cat']  = pd.cut(train['dropoff_latitude'],  n_cat)\n\npickup_masker = pd.pivot_table(train, 'fare_amount', 'pickup_latitude_cat', 'pickup_longitude_cat', 'count')>20\npickup_medians = (pd.pivot_table(train, 'fare_amount', 'pickup_latitude_cat', 'pickup_longitude_cat', 'median')*pickup_masker).fillna(0).sort_index(ascending = False)\ndropoff_masker = pd.pivot_table(train, 'fare_amount', 'dropoff_latitude_cat', 'dropoff_longitude_cat', 'count')>20\ndropoff_medians = (pd.pivot_table(train, 'fare_amount', 'dropoff_latitude_cat', 'dropoff_longitude_cat', 'median')*dropoff_masker).fillna(0).sort_index(ascending = False)\n\nf, axes = plt.subplots(1, 2, figsize = (35, 12))\npickup_hmap = sns.heatmap(pickup_medians, alpha = 0.6, zorder = 2, ax = axes[0], cmap = 'Blues')\npickup_hmap.imshow(nyc_img, aspect = pickup_hmap.get_aspect(),\n                   extent = pickup_hmap.get_xlim() + pickup_hmap.get_ylim(), zorder = 1) \n\npickup_hmap.set_title('Median of fares for pickups[#passengers > 20]', fontsize = 30)\npickup_hmap.set_xlabel('')\npickup_hmap.set_ylabel('')\n\ndropoff_hmap = sns.heatmap(dropoff_medians, alpha = 0.6, zorder = 2, ax = axes[1], cmap = 'Blues')\ndropoff_hmap.imshow(nyc_img, aspect = dropoff_hmap.get_aspect(),\n                    extent = dropoff_hmap.get_xlim() + dropoff_hmap.get_ylim(), zorder = 1) \n\ndropoff_hmap.set_title('Median fares for dropoffs[#passengers > 20]', fontsize = 30)\ndropoff_hmap.set_xlabel('')\ndropoff_hmap.set_ylabel('');","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:06:57.033825Z","iopub.execute_input":"2022-08-10T17:06:57.034340Z","iopub.status.idle":"2022-08-10T17:07:02.849095Z","shell.execute_reply.started":"2022-08-10T17:06:57.034305Z","shell.execute_reply":"2022-08-10T17:07:02.847400Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* This plot strengthens the doubt that travellers picked up or dropped off far from the City have higher fares.\n* One way to capture this is to make new features that will measure the distance and the azimuth from the center to the pickup and dropoff locations (4 features in total).\n* Furthermore adding the interaction terms (azimuth*distance) between them might capture the effect that for example going 10km to the north might not have the same cost as going 10 km to the west. Not sure how great idea it is but the feature importance plot will show us (Spoiler: It's great idea!).\n\n","metadata":{}},{"cell_type":"code","source":"NYC_long = -74.0060\nNYC_lat  = 40.7128\n\ntrain['center_to_pickup_azimuth'], cache, train['center_to_pickup_km'] = geodesic.inv(\n    np.full(train.shape[0], NYC_long), np.full(train.shape[0], NYC_lat),\n    train['pickup_longitude'].values, train['pickup_latitude'].values\n)\ntrain['center_to_dropoff_azimuth'], cache, train['center_to_dropoff_km'] = geodesic.inv(\n    np.full(train.shape[0], NYC_long), np.full(train.shape[0], NYC_lat),\n    train['dropoff_longitude'].values, train['dropoff_latitude'].values\n)\ntest['center_to_pickup_azimuth'], cache, test['center_to_pickup_km'] = geodesic.inv(\n    np.full(test.shape[0], NYC_long), np.full(test.shape[0], NYC_lat),\n    test['pickup_longitude'].values, test['pickup_latitude'].values\n)\ntest['center_to_dropoff_azimuth'], cache, test['center_to_dropoff_km'] = geodesic.inv(\n    np.full(test.shape[0], NYC_long), np.full(test.shape[0], NYC_lat),\n    test['dropoff_longitude'].values, test['dropoff_latitude'].values\n)\n\n#From meters to kilometers\ntrain[['center_to_pickup_km', 'center_to_dropoff_km']] = train[['center_to_pickup_km', 'center_to_dropoff_km']]/1000 \ntest[['center_to_pickup_km', 'center_to_dropoff_km']] = test[['center_to_pickup_km', 'center_to_dropoff_km']]/1000","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:07:13.420620Z","iopub.execute_input":"2022-08-10T17:07:13.421027Z","iopub.status.idle":"2022-08-10T17:07:25.237452Z","shell.execute_reply.started":"2022-08-10T17:07:13.420992Z","shell.execute_reply":"2022-08-10T17:07:25.236385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Datetime features","metadata":{}},{"cell_type":"code","source":"def add_datetime_features(data):\n    data['pickup_datetime'] = pd.to_datetime(data['pickup_datetime'])\n    data['date']            = data['pickup_datetime'].dt.date\n    data['year']            = data['pickup_datetime'].dt.year\n    data['month']           = data['pickup_datetime'].dt.month_name()\n    data['week']            = data['pickup_datetime'].dt.isocalendar().week\n    data['day']             = data['pickup_datetime'].dt.day\n    data['hour']            = data['pickup_datetime'].dt.hour\n    data['day_of_week']     = data['pickup_datetime'].dt.day_name()\n    \n    return data\n\ntrain = add_datetime_features(train).sort_values('pickup_datetime') #Sorting for future plots.\ntest = add_datetime_features(test)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:07:25.240447Z","iopub.execute_input":"2022-08-10T17:07:25.241367Z","iopub.status.idle":"2022-08-10T17:20:21.124428Z","shell.execute_reply.started":"2022-08-10T17:07:25.241317Z","shell.execute_reply":"2022-08-10T17:20:21.122486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.set(rc = {'figure.figsize':(16,6)})\nplt.style.use('seaborn-whitegrid')\nsns.pointplot(x = 'week', y = 'fare_amount', data = train, ci = 95, hue = 'year')\nplt.title('Mean fares by weeks', fontsize = 16);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:22:29.956344Z","iopub.execute_input":"2022-08-10T17:22:29.956894Z","iopub.status.idle":"2022-08-10T17:23:42.133587Z","shell.execute_reply.started":"2022-08-10T17:22:29.956853Z","shell.execute_reply":"2022-08-10T17:23:42.132623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* There are visible differences between the years.\n* Turning date seems to be in the 36th week (The start of September) of the 2012 year. \n* In addition to the year I will add \"after_2012sep\" as a dummy variable.\n\nLet's add a dummy variable and then continue exploration.","metadata":{}},{"cell_type":"code","source":"train['after_2012sep'] = (train['pickup_datetime'] >= '2012-09-01')*1\ntest['after_2012sep'] = (test['pickup_datetime'] >= '2012-09-01')*1","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:23:57.165517Z","iopub.execute_input":"2022-08-10T17:23:57.166914Z","iopub.status.idle":"2022-08-10T17:23:57.210341Z","shell.execute_reply.started":"2022-08-10T17:23:57.166862Z","shell.execute_reply":"2022-08-10T17:23:57.209233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"months = ['January', 'February', 'March', 'April', 'May', 'June',\n          'July', 'August', 'September', 'October', 'November', 'December']\n\nfare_month_day_cross = pd.pivot_table(train, 'fare_amount', 'month', 'day', 'median').reindex(months)\nsns.heatmap(fare_month_day_cross, linewidths = 0.5, cmap=\"coolwarm\")\nplt.title('Median fare each day', fontsize = 16);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:23:59.899110Z","iopub.execute_input":"2022-08-10T17:23:59.899765Z","iopub.status.idle":"2022-08-10T17:24:07.307920Z","shell.execute_reply.started":"2022-08-10T17:23:59.899715Z","shell.execute_reply":"2022-08-10T17:24:07.306362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* From the graph, we observe that the median fare differs by month. At the start of the year median price is low, but till the new year, it increases slowly. Exceptions are holiday months - July and August. \n* Furthermore, there can be seen two holidays which are out of pattern - the start of the new year and Christmas.\n* This graph implies that the month might be an important feature for our prediction.","metadata":{}},{"cell_type":"code","source":"sns.pointplot(x = 'day_of_week', y = 'fare_amount', data = train, ci = 95)\nplt.ylabel('')\nplt.title('Mean fares by weekdays', fontsize = 16);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:26:55.182185Z","iopub.execute_input":"2022-08-10T17:26:55.182592Z","iopub.status.idle":"2022-08-10T17:29:21.786405Z","shell.execute_reply.started":"2022-08-10T17:26:55.182562Z","shell.execute_reply":"2022-08-10T17:29:21.785045Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.pointplot(x = 'day', y = 'fare_amount', data = train, ci = 95)\nplt.ylabel('')\nplt.title('Mean fares by days of months', fontsize = 16);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:29:34.473105Z","iopub.execute_input":"2022-08-10T17:29:34.473717Z","iopub.status.idle":"2022-08-10T17:30:50.789559Z","shell.execute_reply.started":"2022-08-10T17:29:34.473667Z","shell.execute_reply":"2022-08-10T17:30:50.788420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.pointplot(x = 'hour', y = 'fare_amount', data = train, ci = 95, hue = 'day_of_week')\nplt.title('Mean fares by days of months', fontsize = 16);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:32:14.358682Z","iopub.execute_input":"2022-08-10T17:32:14.359050Z","iopub.status.idle":"2022-08-10T17:33:33.718269Z","shell.execute_reply.started":"2022-08-10T17:32:14.359019Z","shell.execute_reply":"2022-08-10T17:33:33.717013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* day doesn't seem a promising feature.\n* hour and weekdays show some patterns.","metadata":{}},{"cell_type":"markdown","source":"# Weather features from external source","metadata":{}},{"cell_type":"markdown","source":"I found weather data for New York on Github. Thanks [zonination](https://github.com/zonination) for sharing it.\n","metadata":{}},{"cell_type":"code","source":"url = 'https://raw.githubusercontent.com/zonination/weather-us/master/nyc.csv'\nweather = pd.read_csv(url, index_col=0)\nweather.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:34:47.017955Z","iopub.execute_input":"2022-08-10T17:34:47.018385Z","iopub.status.idle":"2022-08-10T17:34:47.536791Z","shell.execute_reply.started":"2022-08-10T17:34:47.018349Z","shell.execute_reply":"2022-08-10T17:34:47.535378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* There will be used only those variables which intuitively could have had an effect on fare. ","metadata":{}},{"cell_type":"code","source":"weather['date'] = pd.to_datetime(weather['Date'])\nweather = weather.loc[weather['date'] > '2008-12-31',['date','Mean.TemperatureF','Mean.Wind.SpeedMPH','Events']]\nweather['Events'] = weather['Events'].fillna('unknown')\nweather['Mean.Wind.SpeedMPH'] = weather['Mean.Wind.SpeedMPH'].fillna(weather['Mean.Wind.SpeedMPH'].mean())\nweather['Mean.TemperatureC'] = weather['Mean.TemperatureF'].apply(lambda F: (F-32)*5/9)\nweather['Mean.TemperatureC'] = weather['Mean.TemperatureC'].fillna(weather['Mean.TemperatureC'].mean())\n\nweather = weather.drop('Mean.TemperatureF', axis = 1)\nweather.head()","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:34:49.885209Z","iopub.execute_input":"2022-08-10T17:34:49.885619Z","iopub.status.idle":"2022-08-10T17:34:49.917409Z","shell.execute_reply.started":"2022-08-10T17:34:49.885585Z","shell.execute_reply":"2022-08-10T17:34:49.916090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def merge_weather(data):\n    data['date'] = pd.to_datetime(data['date'])\n    data = data.merge(weather, 'left', 'date')\n    \n    return data\n\ntrain = merge_weather(train)\ntest = merge_weather(test)","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:34:58.768144Z","iopub.execute_input":"2022-08-10T17:34:58.768564Z","iopub.status.idle":"2022-08-10T17:35:12.333254Z","shell.execute_reply.started":"2022-08-10T17:34:58.768506Z","shell.execute_reply":"2022-08-10T17:35:12.331824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(ncols=2, figsize = (15,5))\n\naxs[0].scatter(x = train['Mean.TemperatureC'], y = train['fare_amount'], alpha = 0.1, s = 1)\naxs[0].set_xlabel('TemperatureC')\naxs[0].set_ylabel('Fare')\naxs[0].set_title('Scatterplot of Temperature vs Fare', fontsize = 13)\n\naxs[1].scatter(x = train['Mean.Wind.SpeedMPH'], y = train['fare_amount'], alpha = 0.1, s = 1)\naxs[1].set_xlabel('Wind speed')\naxs[1].set_ylabel('Fare')\naxs[1].set_title('Scatterplot of Wind Speed vs Fare', fontsize = 13);","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2022-08-10T17:42:18.946389Z","iopub.execute_input":"2022-08-10T17:42:18.946842Z","iopub.status.idle":"2022-08-10T17:42:24.076274Z","shell.execute_reply.started":"2022-08-10T17:42:18.946805Z","shell.execute_reply":"2022-08-10T17:42:24.075019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* It seems they both have quadratic relation to fare. ","metadata":{}},{"cell_type":"code","source":"plt.style.use('seaborn-whitegrid')\nfig, axes = plt.subplots(nrows=1,ncols=2, figsize = (23,7))\ntrain['Events'].value_counts(normalize = 'columns').sort_values().plot(kind = 'barh', ax = axes[0])\ntest['Events'].value_counts(normalize = 'columns').sort_values().plot(kind = 'barh', ax = axes[1])\naxes[0].set_title('Distribution of Events [Train]', fontsize = 20)\naxes[1].set_title('Distribution of Events [Test]', fontsize = 20);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:43:12.397803Z","iopub.execute_input":"2022-08-10T17:43:12.398208Z","iopub.status.idle":"2022-08-10T17:43:13.174704Z","shell.execute_reply.started":"2022-08-10T17:43:12.398177Z","shell.execute_reply":"2022-08-10T17:43:13.173648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ax = sns.pointplot(x = 'Events', y = 'fare_amount', data = train, ci = 95)\nax.set_xticklabels(ax.get_xticklabels(), rotation=40, ha=\"right\")\nax.set_title('Mean of fares amount by Events', fontsize = 16);","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:38:02.173761Z","iopub.execute_input":"2022-08-10T17:38:02.174265Z","iopub.status.idle":"2022-08-10T17:40:25.455469Z","shell.execute_reply.started":"2022-08-10T17:38:02.174227Z","shell.execute_reply":"2022-08-10T17:40:25.453807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* One might say that Fog-Thunderstorm and some others are important predictors, but as a portion of them in the dataset is approximately zero, I will omit them.","metadata":{}},{"cell_type":"markdown","source":"# Data preparation","metadata":{}},{"cell_type":"code","source":"train_, val_ = train_test_split(train, train_size = 0.9, random_state = 7)\n\nordinals = ['month','day_of_week'] #They aren't actually ordinal variables but XGBoost doesn't care much.\ninteractions1  = ['travelled_km', 'azimuth']\ninteractions2 = ['pickup_longitude', 'pickup_latitude']\ninteractions3 = ['dropoff_longitude', 'dropoff_latitude']\ninteractions4 = ['center_to_pickup_km', 'center_to_pickup_azimuth']\ninteractions5 = ['center_to_dropoff_km', 'center_to_dropoff_azimuth']\nunchanged = ['year', 'hour', 'Mean.Wind.SpeedMPH', 'Mean.TemperatureC', 'after_2012sep']","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:43:28.968198Z","iopub.execute_input":"2022-08-10T17:43:28.968644Z","iopub.status.idle":"2022-08-10T17:43:39.864719Z","shell.execute_reply.started":"2022-08-10T17:43:28.968608Z","shell.execute_reply":"2022-08-10T17:43:39.863037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"preprocessor = ColumnTransformer(\n    transformers=[\n        ('ordinal', OrdinalEncoder(), ordinals),\n        ('inter1', PolynomialFeatures(interaction_only=True, include_bias = False), interactions1),\n        ('inter2', PolynomialFeatures(interaction_only=True, include_bias = False), interactions2),\n        ('inter3', PolynomialFeatures(interaction_only=True, include_bias = False), interactions3),\n        ('inter4', PolynomialFeatures(interaction_only=True, include_bias = False), interactions4),\n        ('inter5', PolynomialFeatures(interaction_only=True, include_bias = False), interactions5),\n        ('unchanged', 'passthrough', unchanged)\n    ])\n\nX_train = preprocessor.fit_transform(train_)\nX_val = preprocessor.fit_transform(val_)\nX_test = preprocessor.fit_transform(test)\n\ny_train = train_['fare_amount']\ny_val = val_['fare_amount']\n\nX_colnames = \\\nordinals + \\\ninteractions1 + ['(travelled_km*azimuth)'] + \\\ninteractions2 + ['pickup_(longitude*latitude)'] + \\\ninteractions3 + ['dropoff_(longitude*latitude)'] + \\\ninteractions4 + ['center_to_pickup_(km*azimuth)'] + \\\ninteractions5 + ['center_to_dropoff_(km*azimuth)'] + \\\nunchanged","metadata":{"execution":{"iopub.status.busy":"2022-08-10T17:43:42.509734Z","iopub.execute_input":"2022-08-10T17:43:42.510166Z","iopub.status.idle":"2022-08-10T17:43:57.508252Z","shell.execute_reply.started":"2022-08-10T17:43:42.510132Z","shell.execute_reply":"2022-08-10T17:43:57.506853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"__Little comment about interaction1 term:__ In linear regression cutting azimuth variable into categories, one hot encoding and then adding interactions between each one hot encoded variable and distance makes perfect sense. This way you can capture the effect that people going toward the north have a different slope and intercept. For this, you need a lot of dummy variables which is too cumbersome to train. So, instead of doing that, I straightforwardly added interaction between continuous variables. It's not the best way and is hard to interpret even in linear regression but I am planning to use XGBoost and as we will see it does its job very well.\n\n__Same logic applies to another interaction terms.__","metadata":{"execution":{"iopub.status.busy":"2022-08-10T16:10:26.895800Z","iopub.execute_input":"2022-08-10T16:10:26.896377Z","iopub.status.idle":"2022-08-10T16:10:26.933700Z","shell.execute_reply.started":"2022-08-10T16:10:26.896250Z","shell.execute_reply":"2022-08-10T16:10:26.931951Z"}}},{"cell_type":"markdown","source":"# Hyperparameter tuning","metadata":{}},{"cell_type":"markdown","source":"Hyperparameter tuning takes too long. So, I will just provide the best parameters without executing the code.","metadata":{}},{"cell_type":"code","source":"# %%time\n# param_grid = {\n#     'max_depth': [15, 20, 30],\n#     'colsample_bytree' : [0.5, 0.7, 1],\n#     'min_child_weight': [100, 300, 500],\n#     'learning_rate': [0.02, 0.04, 0.1]\n# }\n\n# xgb_gs = GridSearchCV(\n#     xgb.XGBRegressor(n_estimators = 2000, early_stopping_rounds = 30, importance_type = 'gain'),\n#     param_grid,\n#     cv = 2,\n#     scoring = 'neg_root_mean_squared_error',\n#     verbose = 2\n# )\n# xgb_gs.fit(X_train, y_train, eval_set = [(X_val, y_val)])\n# print(xgb_gs.best_params_)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb_mod = xgb.XGBRegressor(max_depth = 30, colsample_bytree = 0.5, learning_rate = 0.04,\n                           min_child_weight = 300, n_estimators = 2000, early_stopping_rounds = 30,\n                           importance_type = 'gain')\nxgb_mod.fit(X_train, y_train, eval_set = [(X_val, y_val)], verbose = False)","metadata":{"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(xgb_mod.feature_importances_, index = X_colnames).sort_values(0)\\\n.plot(kind = 'barh', figsize=(12,6), legend = False)\nplt.title('Feature importances', fontsize = 16);","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"final_preds = xgb_mod.predict(X_test)\nsubmission = pd.DataFrame({'key':test['key'],'fare_amount':final_preds})\nsubmission.to_csv('submission.csv', index = False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"__Suggestions will be appreciated.__","metadata":{"execution":{"iopub.status.busy":"2022-08-10T18:13:34.446419Z","iopub.execute_input":"2022-08-10T18:13:34.446906Z","iopub.status.idle":"2022-08-10T18:13:34.455038Z","shell.execute_reply.started":"2022-08-10T18:13:34.446855Z","shell.execute_reply":"2022-08-10T18:13:34.453758Z"}}}]}