<jupyter_start><jupyter_text>Libraries<jupyter_code>import folium
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from glmnet import ElasticNet
from scipy.stats import zscore
from sklearn.utils import resample

import glob
import json
import os
import warnings
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>----------<jupyter_code>%time taxi_df = pd.read_feather("../preprocessed_data/large/yellow_tripdata.feather")
taxi_df.dtypes
taxi_df.tail()<jupyter_output><empty_output><jupyter_text>----------<jupyter_code># binning methods (from lab-2)
def sturges(x):
    from numpy import log2
    return int(log2(x)) + 1

def rice(x):
    return int(2 * x ** (1/3))

def scott(large, small, sd, x):
    return int((large - small) / (3.5 * (sd/x ** (1/3))))

def freedman(large, small, iqr, x):
    return int((large - small) / (2 * (iqr/x ** (1/3))))

def square(x):
    from numpy import sqrt
    return int(sqrt(x))
def distplot_all_bins(col):
    """
        col is the series of data to plot
    """
    MAX, MIN, SD, IQR, N = col.max(), col.min(), col.std(), col.quantile(0.75)-col.quantile(0.25), len(col)
    bins_by_title = {"Sturges Binnings": sturges(N),
                     "Rice Binnings": rice(N),
                     "Scott Binnings": scott(MAX,MIN,SD,N),
                     "Freedman Binnings": freedman(MAX,MIN,IQR,N),
                     "Square Binnings": square(N)}
    figs = []    
    for t in bins_by_title:
        print(bins_by_title[t], "bins")
        %time figs.append(sns.distplot(col, bins=bins_by_title[t]))
        plt.title(t)
        plt.show()
    
    return figs
def logify(x):
    from numpy import log
    return log(x) if x else 0<jupyter_output><empty_output><jupyter_text>----------  Map Visualisation<jupyter_code># merge with taxi zone geometry
sf = gpd.read_file("../raw_data/taxi_zones.shp").drop(["OBJECTID","Shape_Leng","zone","borough"], axis=1)
# convert the geometry shape to latitude and longitude (from lab-2)
sf["geometry"] = sf["geometry"].to_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")<jupyter_output><empty_output><jupyter_text>A driver cannot decide the destination of a trip, so on behalf of a driver, we should focus on pickup location rather than dropoff.<jupyter_code>zone_mean_tip_ratio_df = taxi_df[["PULocationID","tip_ratio"]].groupby("PULocationID").mean().reset_index()
zone_mean_tip_ratio_df = gpd.GeoDataFrame(pd.merge(zone_mean_tip_ratio_df, sf, left_on="PULocationID", right_on="LocationID")).drop("LocationID", axis=1)
geoJSON = zone_mean_tip_ratio_df[["PULocationID","geometry"]].to_json()<jupyter_output><empty_output><jupyter_text>Distributional plots.<jupyter_code>figs = distplot_all_bins(zone_mean_tip_ratio_df["tip_ratio"])
figs[4].get_figure().savefig("../plots/DistplotMeanTipRatio.png")<jupyter_output><empty_output><jupyter_text>**Plot tip_ratio.**<jupyter_code>m_tip_ratio = folium.Map(location=[40.66,-73.94], tiles="Stamen Terrain", zoom_start=10)

folium.Choropleth(geo_data=geoJSON,
                  name="choropleth",
                  data=zone_mean_tip_ratio_df,
                  columns=["PULocationID", "tip_ratio"],
                  key_on="properties.PULocationID",
                  fill_color="OrRd",
                  fill_opacity=0.9,
                  line_opacity=0.5,
                  legend_name="Tip Ratio"
                 ).add_to(m_tip_ratio)

m_tip_ratio.save("../plots/FoliumChoroplethMapMeanTipRatio.html")
m_tip_ratio<jupyter_output><empty_output><jupyter_text>---------- Apart from high tip ratio, a driver's income will also increase if there are- large fare amount each trip, and- many trips in a small area of region.<jupyter_code>zone_sum_fare_df = taxi_df[["PULocationID","fare_amount"]].groupby("PULocationID").sum().reset_index()
zone_sum_fare_df = gpd.GeoDataFrame(pd.merge(zone_sum_fare_df, sf, left_on="PULocationID", right_on="LocationID")).drop("LocationID", axis=1)
geoJSON = zone_sum_fare_df[["PULocationID","geometry"]].to_json()

zone_sum_fare_df["fare_per_area"] = zone_sum_fare_df["fare_amount"] / zone_sum_fare_df["Shape_Area"]<jupyter_output><empty_output><jupyter_text>Distributional plots.<jupyter_code>figs = distplot_all_bins(zone_sum_fare_df["fare_per_area"])
figs[0].get_figure().savefig("../plots/DistplotSumFarePerArea.png")<jupyter_output><empty_output><jupyter_text>Perform log-transformation.<jupyter_code>zone_sum_fare_df["fare_per_area_log"] = zone_sum_fare_df["fare_per_area"].apply(logify)
col = zone_sum_fare_df["fare_per_area_log"]
b = sturges(len(col))
print(b, "bins")
sns.distplot(col, bins=b)
plt.title("Sturges Binnings")
plt.savefig("../plots/DistplotSumFarePerAreaLog.png")
plt.show()<jupyter_output>9 bins<jupyter_text>**Plot log(fare_per_area).**<jupyter_code>m_fare_per_area_log = folium.Map(location=[40.66,-73.94], tiles="Stamen Terrain", zoom_start=10)

folium.Choropleth(geo_data=geoJSON,
                  name="choropleth",
                  data=zone_sum_fare_df,
                  columns=["PULocationID", "fare_per_area_log"],
                  key_on="properties.PULocationID",
                  fill_color="OrRd",
                  fill_opacity=0.9,
                  line_opacity=0.5,
                  legend_name="Log of Fare Per Area"
                 ).add_to(m_fare_per_area_log)

m_fare_per_area_log.save("../plots/FoliumChoroplethMapSumFarePerAreaLog.html")
m_fare_per_area_log<jupyter_output><empty_output><jupyter_text>----------  Distributional Plots of Features Fare amount:<jupyter_code>figs = distplot_all_bins(taxi_df["fare_amount"])
figs[0].get_figure().savefig("../plots/DistplotFareAmount.png")<jupyter_output><empty_output><jupyter_text>Tolls amount:<jupyter_code>figs = distplot_all_bins(taxi_df.loc[taxi_df["tolls_amount"] != 0, "tolls_amount"])
figs[0].get_figure().savefig("../plots/DistplotTollsAmount.png")<jupyter_output><empty_output><jupyter_text>Trip distance:<jupyter_code>figs = distplot_all_bins(taxi_df["trip_distance"])
figs[1].get_figure().savefig("../plots/DistplotTripDistance.png")<jupyter_output><empty_output><jupyter_text>Trip time:<jupyter_code>figs = distplot_all_bins(taxi_df["trip_time"])
figs[2].get_figure().savefig("../plots/DistplotTripTime.png")<jupyter_output><empty_output><jupyter_text>Average speed:<jupyter_code>figs = distplot_all_bins(taxi_df["avg_speed"])
figs[1].get_figure().savefig("../plots/DistplotAvgSpeed.png")<jupyter_output><empty_output><jupyter_text>Tip ratio:<jupyter_code>figs = distplot_all_bins(taxi_df["tip_ratio"])
figs[0].get_figure().savefig("../plots/DistplotTipRatio.png")<jupyter_output><empty_output><jupyter_text>----------  Feature Selection Assume linear relationships between features, then:<jupyter_code>%%time
CORR_COLS = ["passenger_count", "RatecodeID", "trip_distance", "trip_time", "avg_speed", "fare_amount", \
             "tolls_amount", "rushhour", "overnight", "weekend_or_holiday", "PRCP", "SNOW", "TMAX", "TMIN", "tip_ratio"]
fig, ax = plt.subplots(figsize=(15,10))
sns.heatmap(taxi_df[CORR_COLS].corr(method="pearson"), annot=True, vmin=-1, vmax=1, center=0)
plt.savefig("../plots/HeatmapNumericFeatures.png")
plt.show()<jupyter_output><empty_output><jupyter_text>Perform log-transformation on right-skewed features.<jupyter_code>%%time
taxi_df["fare_amount_log"] = taxi_df["fare_amount"].apply(logify)
taxi_df["trip_distance_log"] = taxi_df["trip_distance"].apply(logify)
taxi_df["trip_time_log"] = taxi_df["trip_time"].apply(logify)
taxi_df["avg_speed_log"] = taxi_df["avg_speed"].apply(logify)
%%time
CORR_COLS = ["passenger_count", "RatecodeID", "trip_distance_log", "trip_time_log", "avg_speed_log", "fare_amount_log", \
             "tolls_amount", "rushhour", "overnight", "weekend_or_holiday", "PRCP", "SNOW", "TMAX", "TMIN", "tip_ratio"]
fig, ax = plt.subplots(figsize=(15,10))
sns.heatmap(taxi_df[CORR_COLS].corr(method="pearson"), annot=True, vmin=-1, vmax=1, center=0)
plt.savefig("../plots/HeatmapNumericFeaturesLog.png")
plt.show()<jupyter_output><empty_output><jupyter_text>----------<jupyter_code>%%time
taxi_df[CORR_COLS[:-1]] = taxi_df[CORR_COLS[:-1]].apply(zscore)   # standardisation
X = taxi_df[CORR_COLS[:-1]]
y = taxi_df[CORR_COLS[-1]]

lasso_fit = ElasticNet()
lasso_fit.fit(X.values, y.values)
print("Best lambda value for LASSO:", {lasso_fit.lambda_best_[0]})

pd.DataFrame(index = ["Intercept"] + CORR_COLS[:-1],
             data = [lasso_fit.intercept_] + list(lasso_fit.coef_), columns = ["Coefficient"])<jupyter_output>Best lambda value for LASSO: {7.45905961901867e-05}
CPU times: user 4min 50s, sys: 4min 44s, total: 9min 35s
Wall time: 11min 51s<jupyter_text>This makes the model to exclude features: passenger_count, avg_speed_log, overnight, PRCP and TMAX.<jupyter_code>selected = list(set(CORR_COLS[:-1])-set(["passenger_count","avg_speed_log","overnight","PRCP","TMAX"]))
selected<jupyter_output><empty_output><jupyter_text>----------  Binning up tip_ratio The suggested tipping rate for taxi drivers in NYC is 10%-20% (of the bill).https://www.fodors.com/world/north-america/usa/new-york/new-york-city/travel-tips/tipping-52564530Therefore, we split tip_ratio into 3 ranges, low (20%).<jupyter_code>%%time
model_df = taxi_df[selected]
model_df["_tip_ratio"] = taxi_df["tip_ratio"]
model_df.loc[model_df["_tip_ratio"] < 0.1, "tip_ratio"] = "l"   # low
model_df.loc[(model_df["_tip_ratio"] >= 0.1) & (model_df["_tip_ratio"] <= 0.2), "tip_ratio"] = "m"   # moderate
model_df.loc[model_df["_tip_ratio"] > 0.2, "tip_ratio"] = "h"   # high
model_df = model_df.drop("_tip_ratio", axis=1)
model_df.tail()
model_df["tip_ratio"].value_counts().to_dict()
model_df.describe()<jupyter_output><empty_output><jupyter_text>----------  Resampling 10% of the original size, in a stratified fashion.<jupyter_code>%time model_df = resample(model_df, replace=True, n_samples=len(model_df)/10, random_state=30034, stratify=model_df["tip_ratio"])
model_df<jupyter_output>CPU times: user 21 s, sys: 524 ms, total: 21.5 s
Wall time: 21.4 s<jupyter_text>----------  Save to `.feather`<jupyter_code>%time model_df.reset_index().drop("index", axis=1).to_feather("../preprocessed_data/large/yellow_tripmodel.feather")<jupyter_output>CPU times: user 584 ms, sys: 133 ms, total: 717 ms
Wall time: 439 ms
{'fraction_non_alphanumeric': 0.15868699621248908, 'fraction_numerical': 0.01699524133242692, 'mean_word_length': 5.432229856339788, 'pattern_counts': {'":': 5, '<': 95, '<?xml version=': 0, '>': 95, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '518424', 'license': 'MIT', 'n_tokens_mistral': 4060, 'n_tokens_neox': 3843, 'n_words': 577, 'path': 'code/3_Analysis.ipynb', 'repo_name': 'zzhcai/NYCTaxi', 'size': '10297'}