<jupyter_start><jupyter_text>New York City Taxi Analysis--- IntroductionThis notebook shows you how we can use exploratory data analysis and machine learning for predicting breast cancers. All our models are built using  New York City Taxi dataset.Also, based on the data in this dataset, we would like to suggest following projects which create valuable business opportunities for **NYC Taxi and Limousine Commission** (**TLC**).1. Intelligent Cab Positioning: position cabs where they are most needed, how many taxis to dispatch, and how ridership varies over time. 2.  Predicting the Number of Taxi Pickups: given a one-hour time window and a location within New York City3. Real-time Allocation Suggestion for Drivers: Based on the past data available to us, we can create a model which can suggest most profitable areas for cab drivers to pick customers.4. Predicting Travel Time for Each Taxi Ride: This helps TLC to dispatch cabs to different locations in the city. In addition to that this information helps the customer to plan their travels as well.In this node book, we are going to analyze item number four (i.e. **Predicting Travel Time for Each Taxi Ride**)  Loading Necessary Tools (Python/Numpy, Pandas, Scikit-Learn, LGBM, and etc.)<jupyter_code>import numpy as np
import pandas as pd
import lightgbm as lgb
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import KFold
from sklearn.metrics import confusion_matrix
from sklearn.preprocessing import scale

sns.set()

print('Numpy version: {}'.format(np.__version__))
print('Pandas version: {}'.format(pd.__version__))
print('Light Gradient Booting version: {}'.format(lgb.__version__))

%matplotlib inline

import warnings
warnings.filterwarnings('ignore')<jupyter_output>Numpy version: 1.13.3
Pandas version: 0.20.3
Light Gradient Booting version: 2.1.0<jupyter_text>Loading the DatasetNext, we are going to load the dataset as a Pandas data frame (DataFrame). In ML model building, it is the best practice to keep the testing dataset aside and it should not be mixed with training dataset. So as the very first step, we split our dataset into training and testing subsets.We are going to use **80 percent** (i.e. **8000 data rows**) of the dataset for training and the rest (i.e. **2000 data rows**) for testing.<jupyter_code>data = pd.read_csv('./data/yellow_tripdata_2016-01_10000rows.csv')
data = data.sample(n=len(data))

train_size = int(len(data)*0.8)
test_size = len(data) - train_size

assert len(data) == (train_size + test_size)
print('training dataset size: {}'.format(train_size))
print('testing dataset size: {}'.format(test_size))

train = data.iloc[0:train_size]
test = data.iloc[train_size:]

assert len(data) == len(train) + len(test)<jupyter_output>training dataset size: 8000
testing dataset size: 2000<jupyter_text>As you know real-world datasets come with missing values  (or measurements). So, first of all, we need to check this dataset contains null or missing values. If it has missing values, proper data imputation techniques need to be applied.I'm going to use **Panda's** **`isnull()`** for testing null/missing values.<jupyter_code>assert train.isnull().values.any() == False
assert test.isnull().values.any() == False<jupyter_output><empty_output><jupyter_text>Great, this dataset doesn't have any missing values so we don't need to use data imputation techniques. If it contains missing values I might be using one or more following techniques.1. **Impute missing values with mean**: Usually, this is suitable for numerical features.2. **Impute missing values with median**: This approach is more appropriate for categorical features.3. **EM algorithm based approaches**: Bit advanced but very useful tool for missing value imputation.  Create Response VariableIn this analysis, we are going to predict taxi trip duration using available data. So we create our response variable using **tpep_pickup_datetime** and **tpep_dropoff_datetime**.<jupyter_code>from datetime import datetime

fmt = '%m/%d/%Y %H:%M'

def calculate_trip_duration(row):
    start = datetime.strptime(row['tpep_pickup_datetime'], fmt)
    end = datetime.strptime(row['tpep_dropoff_datetime'], fmt)
    return (end - start).seconds / 60.0

train['trip_time'] = train.apply(calculate_trip_duration, axis=1)
test['trip_time'] = test.apply(calculate_trip_duration, axis=1)<jupyter_output><empty_output><jupyter_text>Exploratory Data Analysis and VisualizationExploratory data analysis helps us to understand the nature of the training dataset. Also, it is a great tool to understand the relationships between different features of the dataset.  Distribution of  Response Variable (Taxi Trip Duration)As the first step of our exploratory data analysis, we plot the distribution of taxi trip durations measured in minutes.<jupyter_code>f, axes = plt.subplots(1, 1)
f.set_size_inches(7, 5)

sns.boxplot(train['trip_time']).set_title("Distrubution of Trip Time")
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Above graph clearly shows we have few **outliers**. So let's remove those and redraw the graph again.<jupyter_code>train = train.drop(train[train.trip_time > 120].index)
test = test.drop(test[test.trip_time > 120].index)
f, axes = plt.subplots(1, 1)
f.set_size_inches(7, 5)

#sns.boxplot(test['trip_time'])
sns.boxplot(train['trip_time']).set_title("Distrubution of Trip Time")
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Now response variables look OK.  Various Plots to Identify Good Features and Relationship Between FeaturesWe can use the following three plots to identify quality (in terms of predictive power) features for our model building process. Also, it is worth to note that, for real-world projects, we might use a lot more plots and visualizations than we presented in this notebook.<jupyter_code>f, axes = plt.subplots(1, 1)
f.set_size_inches(7, 5)

sns.boxplot(train['trip_distance']).set_title("Distrubution of Trip Distance")
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Scatter Plots: Pickup and Dropoff Points<jupyter_code>f, axes = plt.subplots(1, 2)
f.set_size_inches(12, 7)
axes[0].scatter(train.pickup_longitude, train.pickup_latitude)
axes[0].set(title='Distribution of Pickup Points', 
            xlabel='pickup longitude', ylabel='pickup latitude')

axes[1].scatter(train.dropoff_longitude, train.dropoff_latitude)
axes[1].set(title='Distribution of Dropoff Points', 
            xlabel='dropoff longitude', ylabel='dropoff latitude')

plt.show()<jupyter_output><empty_output><jupyter_text>So it looks like we have few **outliers** based on the longitude and altitude data. So let's further clean our dataset.<jupyter_code>cond_1 = (train.dropoff_longitude > -73.5) 
cond_2 = (train.dropoff_latitude < 40.5 )
train = train.drop(train[cond_1 & cond_2].index)
#test = test.drop(test[cond_1 & cond_2].index)

cond_3 = (train.dropoff_longitude > -73.5) 
cond_4 = (train.dropoff_latitude < 40.5 )
train = train.drop(train[cond_3 & cond_4].index)
#test = test.drop(test[cond_1 & cond_2].index)
f, axes = plt.subplots(1, 2)
f.set_size_inches(12, 7)
axes[0].scatter(train.pickup_longitude, train.pickup_latitude)
axes[0].set(title='Distribution of Pickup Points', 
            xlabel='pickup longitude', ylabel='pickup latitude')

axes[1].scatter(train.dropoff_longitude, train.dropoff_latitude)
axes[1].set(title='Distribution of Dropoff Points', 
            xlabel='dropoff longitude', ylabel='dropoff latitude')

plt.show()<jupyter_output><empty_output><jupyter_text>Great, now data looks OK in terms of longitudes and altitudes   The relationship between Longitude/Latitude vs Trip Time.We use following **3D plots** to check the relationship between our response variable and longitudes/latitudes.<jupyter_code>from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
fig.set_size_inches(13, 4)

ax = fig.add_subplot(121, projection='3d')
ax.scatter(train.pickup_latitude.values,
            train.pickup_longitude.values,
            train.trip_time.values)
ax.set_xlabel('pickup_latitude')
ax.set_ylabel('pickup_longitude')
ax.set_zlabel('trip_time')


ax = fig.add_subplot(122, projection='3d')
ax.scatter(train.dropoff_latitude.values,
            train.dropoff_longitude.values,
            train.trip_time.values)

ax.set_xlabel('dropoff_latitude')
ax.set_ylabel('dropoff_longitude')
ax.set_zlabel('trip_time')

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>OK, it looks like there is **no apparent relationship** between trip tims and longitudes/latitudes. So  I'm going to engineer some features out of longitudes and latitudes.<jupyter_code>AVG_EARTH_RADIUS = 6371  # in km

def haversine(point1, point2, miles=False):
    """ Calculate the great-circle distance between two points on the Earth surface.
    :input: two 2-tuples, containing the latitude and longitude of each point
    in decimal degrees.
    Example: haversine((45.7597, 4.8422), (48.8567, 2.3508))
    :output: Returns the distance bewteen the two points.
    The default unit is kilometers. Miles can be returned
    if the ``miles`` parameter is set to True.
    """

    # Unpack latitude/longitude
    lat1, lng1 = point1
    lat2, lng2 = point2

    # Convert all latitudes/longitudes from decimal degrees to radians
    lat1, lng1, lat2, lng2 = map(np.radians, (lat1, lng1, lat2, lng2))

    # Calculate haversine
    lat, lng= lat2 - lat1, lng2 - lng1
    d = np.sin(lat * 0.5) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(lng * 0.5) ** 2
    h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d))
    return h * 0.621371 if miles else h # in m/km

def dummy_manhattan_distance(point1, point2):
    lat1, lng1 = point1
    lat2, lng2 = point2
    a = haversine((lat1, lng1), (lat1, lng2))
    b = haversine((lat1, lng1), (lat2, lng1))

    return haversine((lat1, lng1), (lat1, lng2)) + haversine((lat1, lng1), (lat2, lng1))
# apply calculated distances to both training and tesing 
# data sets.

# training
train['distance_haversine'] = haversine(
    (train['pickup_latitude'].values, train['pickup_longitude'].values),
    (train['dropoff_latitude'].values, train['dropoff_longitude'].values)
)

train['distance_dummy_manhattan'] = dummy_manhattan_distance(
    (train['pickup_latitude'].values, train['pickup_longitude'].values),
    (train['dropoff_latitude'].values, train['dropoff_longitude'].values)
)

# testing 
test['distance_haversine'] = haversine(
    (test['pickup_latitude'].values, test['pickup_longitude'].values),
    (test['dropoff_latitude'].values, test['dropoff_longitude'].values)
)

test['distance_dummy_manhattan'] = dummy_manhattan_distance(
    (test['pickup_latitude'].values, test['pickup_longitude'].values),
    (test['dropoff_latitude'].values, test['dropoff_longitude'].values)
)<jupyter_output><empty_output><jupyter_text>In the following graph, I'm testing my new features against the predictive variable i.e. **trip time**.<jupyter_code>independent = ['distance_haversine', 'distance_dummy_manhattan']
dependents =  ['trip_time']

f, axes = plt.subplots(1, 2)
f.set_size_inches(12, 7)

for i in range(len(independent)):
    sns.regplot(x=independent[i], 
                y=dependents[0], data=train, ax=axes[i])

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Great, now we have engineered **two features** which can be **really useful** when it comes to **model building**.  Additionally, we are going to create several more features as described below.All these features are based on **`tpep_pickup_datetime`**. So first of all, let's convert **`tpep_pickup_datetime`** to a **`datetime`** series.<jupyter_code>train['pickup_datetime'] = pd.to_datetime(train['tpep_pickup_datetime'])
test['pickup_datetime'] = pd.to_datetime(train['tpep_pickup_datetime'])<jupyter_output><empty_output><jupyter_text>Now, I'm going to create features such as **pickup_weekday**, **pickup_hour**, **night_trip**, and etc.<jupyter_code>#astype('category')
#['pickup_weekday', 'rush_hour', 'weekday', 
               #'night_trip', 'pickup_day', 'payment_type']
#print(type(data['payment_type']))

train['pickup_date'] = train['pickup_datetime'].dt.date

train['pickup_weekday'] = train['pickup_datetime'].dt.weekday
train['pickup_weekday'] = train['pickup_weekday'].astype('category')


train['pickup_day'] = train['pickup_datetime'].dt.day
train['pickup_day'] = train['pickup_day'].astype('category')

train['pickup_month'] = train['pickup_datetime'].dt.month
train['pickup_hour'] = train['pickup_datetime'].dt.hour
train['pickup_minute'] = train['pickup_datetime'].dt.minute
train['pickup_am'] = train['pickup_hour'] < 12

train['night_trip'] = [True if x < 7 else False for x in train['pickup_hour']]
train['night_trip'] = train['night_trip'].astype('category')

train['rush_hour'] = [True if 9 < x < 20 else False for x in train['pickup_hour']]
train['rush_hour'] = train['rush_hour'].astype('category')

train['weekday'] = [True if x < 5 else False for x in train['pickup_weekday']]
train['weekday']  = train['weekday'].astype('category')

train['payment_type'] = train['payment_type'].astype('category')



test['pickup_date'] = test['pickup_datetime'].dt.date

test['pickup_weekday'] = test['pickup_datetime'].dt.weekday
test['pickup_weekday'] = test['pickup_weekday'].astype('category')

test['pickup_day'] = test['pickup_datetime'].dt.day
test['pickup_day'] = test['pickup_day'].astype('category')

test['pickup_month'] = test['pickup_datetime'].dt.month
test['pickup_hour'] = test['pickup_datetime'].dt.hour
test['pickup_minute'] = test['pickup_datetime'].dt.minute
test['pickup_am'] = test['pickup_hour'] < 12

test['night_trip'] = [True if x < 7 else False for x in test['pickup_hour']]
test['night_trip'] = test['night_trip'].astype('category')

test['rush_hour'] = [True if 9 < x < 20 else False for x in test['pickup_hour']]
test['rush_hour'] = test['rush_hour'].astype('category')

test['weekday'] = [True if x < 5 else False for x in test['pickup_weekday']]
test['weekday']  = test['weekday'].astype('category')

test['payment_type'] = test['payment_type'].astype('category')<jupyter_output><empty_output><jupyter_text>Now let's plot some of these **new features** against our **response variable**.<jupyter_code>import calendar

f, axes = plt.subplots(1, 1)
f.set_size_inches(16, 8)
sns.pointplot(x='pickup_hour',y='trip_time',data=train,
             hue='pickup_weekday', ci=None).set_title('Trip Time vs Pickup Hour')

plt.show()
independent = ['rush_hour', 'night_trip', 'weekday']
dependents =  ['trip_time']

f, axes = plt.subplots(1, 3)
f.set_size_inches(10, 4)

for i in range(len(independent)):
    sns.barplot(x=independent[i], y=dependents[0], data=train, ax=axes[i])

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>According to above plots, it is obvious that these generated features might be good as features when it comes to building our models.**Note**: We just plotted a couple of graphs only. But in practice, it is worth to check the relationship between all possible feature combinations.  Machine Learning Model BuildingNow its time to build some machine learning models. We will start with simple models with few features. But as we progress through the model building process we will explore more advanced models as well. Linear RegressionAs the initial step, we will create a Linear Regression model just using **trip_distance** only. So this model will work as our baseline model. Also, it is worth to note that, we are using **`5-Fold`**  cross-validation and **`root mean squared error`** as our evaluation metric.<jupyter_code>%%time
X_train = train[['trip_distance']].values
y_train = train['trip_time'].values

kf = KFold(n_splits=5, shuffle=True, random_state=1024)

cv_accuracies = []
for train_index, test_index in kf.split(X_train):
    X_train_cv = X_train[train_index]
    y_train_cv = y_train[train_index]
    
    X_test_cv = X_train[test_index]
    y_test_cv = y_train[test_index]
    
    lr_simple = LinearRegression()
    lr_simple.fit(scale(X_train_cv), y_train_cv)
    y_predict_cv = lr_simple.predict(scale(X_test_cv))
    current_accuracy = np.sqrt(mean_squared_error(y_test_cv, y_predict_cv))
    cv_accuracies.append(current_accuracy)
    print('root mean squared error: %f' %(current_accuracy))
    
print ('---------------------------------------')
print ('average root mean squared error: %f' %(sum(cv_accuracies)/len(cv_accuracies)))<jupyter_output>root mean squared error: 6.670934
root mean squared error: 6.281705
root mean squared error: 5.865452
root mean squared error: 6.008910
root mean squared error: 6.172982
---------------------------------------
average root mean squared error: 6.199996
CPU times: user 15.6 ms, sys: 0 ns, total: 15.6 ms
Wall time: 18.3 ms<jupyter_text>Well, as the very starting model, the result looks  OK. But, let's add few more features and look the predictive power of our model.So I'm going to use following numerical and categorical features. All categorical variables will be converted to **one-hot** vectors using **`pd.get_dummies`** method.<jupyter_code>cat_columns = ['rush_hour', 'weekday', 'pickup_weekday',
               'night_trip', 'payment_type', 'pickup_day']

numerical_columns = ['trip_distance', 'distance_haversine', 'distance_dummy_manhattan']

train_cat_onehot =  pd.get_dummies(train[cat_columns])
X_train  = pd.concat([train[numerical_columns],
                     train_cat_onehot], axis=1).values

test_cat_onehot =  pd.get_dummies(test[cat_columns])
test_cat_onehot = test_cat_onehot.reindex(columns=train_cat_onehot.columns,
                                        fill_value=0)

X_test  = pd.concat([test[numerical_columns],
                      test_cat_onehot], axis=1).values<jupyter_output><empty_output><jupyter_text>OK, we can see some improvement, unfortunately, that is not significant. So next, I'm going to use few more ML models to check whether I can further improve **`RMSE`**.  Random Forest<jupyter_code>%%time

kf = KFold(n_splits=5, shuffle=True, random_state=1024)

cv_accuracies = []
for train_index, test_index in kf.split(X_train):
    X_train_cv = X_train[train_index]
    y_train_cv = y_train[train_index]
    
    X_test_cv = X_train[test_index]
    y_test_cv = y_train[test_index]
    
    lr_simple = RandomForestRegressor(n_estimators=1000, max_depth=3)
    lr_simple.fit(scale(X_train_cv), y_train_cv)
    y_predict_cv = lr_simple.predict(scale(X_test_cv))
    current_accuracy = np.sqrt(mean_squared_error(y_test_cv, y_predict_cv))
    cv_accuracies.append(current_accuracy)
    print('root mean squared error: %f' %(current_accuracy))
    
print ('---------------------------------------')
print ('average root mean squared error: %f' %(sum(cv_accuracies)/len(cv_accuracies)))<jupyter_output>root mean squared error: 6.138930
root mean squared error: 5.593333
root mean squared error: 5.396727
root mean squared error: 5.448196
root mean squared error: 5.562827
---------------------------------------
average root mean squared error: 5.628003
CPU times: user 31.3 s, sys: 78.1 ms, total: 31.4 s
Wall time: 31.3 s<jupyter_text>Light Gradient Booting Regression<jupyter_code>%%time

kf = KFold(n_splits=5, shuffle=True, random_state=1024)

cv_accuracies = []
for train_index, test_index in kf.split(X_train):
    X_train_cv = X_train[train_index]
    y_train_cv = y_train[train_index]
    
    X_test_cv = X_train[test_index]
    y_test_cv = y_train[test_index]
    
    lr_simple = lgb.LGBMRegressor(n_estimators=500)
    lr_simple.fit(scale(X_train_cv), y_train_cv)
    y_predict_cv = lr_simple.predict(scale(X_test_cv))
    current_accuracy = np.sqrt(mean_squared_error(y_test_cv, y_predict_cv))
    cv_accuracies.append(current_accuracy)
    print('root mean squared error: %f' %(current_accuracy))
    
print ('---------------------------------------')
print ('average root mean squared error: %f' %(sum(cv_accuracies)/len(cv_accuracies)))

%%time

from sklearn.grid_search import GridSearchCV
from sklearn.metrics.scorer import make_scorer

def MSE(y_true,y_pred):
    rmse = np.sqrt(mean_squared_error(y_true, y_pred))
    return rmse

my_scorer = make_scorer(MSE, greater_is_better=False)

#cls = lgb.LGBMRegressor()
cls = RandomForestRegressor()

param_grid = {
    'max_depth': [7, 13,],
    'n_estimators': [50, 100, 150, 500]
}

gs = GridSearchCV(estimator=cls, param_grid=param_grid, 
                  cv=5, n_jobs=8, scoring=my_scorer, verbose = 0)
gs.fit(X_train, y_train)<jupyter_output>CPU times: user 15.6 s, sys: 531 ms, total: 16.2 s
Wall time: 1min 23s<jupyter_text>OK, now its time to do some hyper-parameter optimization. When it comes to the hyper-parameter optimization we can use:1. Grid Search2. Random Search 3. Baysean Hyper-Parameter OptimizationFor this project, we are going to use **`Grid Search`**  implementation which comes with **`Scikit-Learn`**.<jupyter_code>%%time

kf = KFold(n_splits=5, shuffle=True, random_state=1024)

cv_accuracies = []
for train_index, test_index in kf.split(X_train):
    X_train_cv = X_train[train_index]
    y_train_cv = y_train[train_index]
    
    X_test_cv = X_train[test_index]
    y_test_cv = y_train[test_index]
    
    lr_simple = gs.best_estimator_
    lr_simple.fit(scale(X_train_cv), y_train_cv)
    y_predict_cv = lr_simple.predict(scale(X_test_cv))
    current_accuracy = np.sqrt(mean_squared_error(y_test_cv, y_predict_cv))
    cv_accuracies.append(current_accuracy)
    print('root mean squared error: %f' %(current_accuracy))
    
print ('---------------------------------------')
print ('average root mean squared error: %f' %(sum(cv_accuracies)/len(cv_accuracies)))<jupyter_output>root mean squared error: 5.802245
root mean squared error: 5.216398
root mean squared error: 5.042038
root mean squared error: 4.961596
root mean squared error: 5.224501
---------------------------------------
average root mean squared error: 5.249356
CPU times: user 56.6 s, sys: 78.1 ms, total: 56.7 s
Wall time: 57.1 s<jupyter_text>Finally, we are ready to run our best model against the **testing** dataset.We run our best model against the test dataset and report the **RMSE** value.<jupyter_code>y_test = test['trip_time'].values
y_pred = gs.best_estimator_.predict(scale(X_test))
print('Root mean squared error of training dataset: {0:.3f} '.format(np.sqrt(mean_squared_error(y_test, y_pred))))<jupyter_output>Root mean squared error of training dataset: 6.091<jupyter_text>Finally, we plot sample predictions against actual trip time as shown below.<jupyter_code>f, (ax1, ax2) = plt.subplots(2, 1, sharey=True)
f.set_size_inches(12, 4)

s1 = y_test[700:750]
ax1.bar(np.arange(len(s1)), s1)
ax1.set_title('Actual Trip Time')

s1 = y_pred[700:750]
ax2.bar(np.arange(len(s1)), s1, color='r')
ax2.set_title('Predicted Trip Time')

plt.tight_layout()
plt.show()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup