<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Use the bash =)
Step1:  Rename the columns. They're now streamlined with the taxi input.
Step2:  Parse the dates.
Step3:  Have a look at every station. How many pickups are there per station?
Step4:  When you take a look at the map, you won't be surprised that the Grand Central Station is the one with the most frequent usage.
Step5:  So we have 326 different start stations.
Step6:  We have 326 different end stations, too. But are all end_stations contained in start_stations and vive versa?
Step7:  So this is the proof, that all stations are used as start and end stations.
Check for missing and false data
Step8:  So there is not that much data missing. That's quite surprising, maybe it's wrong.
Step9:  So we only have many zeros in the gender-feature. Gender can have 3 values (0=unknown, 1=male, 2=female)
Step10:  So we see, we have also missing values in "gender" and "birth year". But we are not interested in this features, therefore we can ignore this missing values at this point.
Quick pverview about the trip_times
Step11:  We are interested in the trip time in minutes.
Step12:  We have lots of outliers.
So as we can see, we have many outliers.
Step13:  Identify the the cases without geo data and remove them from our data to be processed.
Step14:  So how many percent of data are left to be processed?
Step15:  <font color = 'blue'> So we dropped nothing because of missing geo tags. </font color>
Step16:  We sometimes have some unreasonably small trip_times.
Calculate the distance of a trip. Use the difference of lon / lat data and parse it to metrics.
Step17:  Be aware, this operation is quite slow.
Step18:  We cannot use miles because vincenty ceils / floors the results
Step19:  Check if new column is present
Step20:  Parse to km/h
Step21:  Convert avg_velocity from km/h to miles/h. (1km = 0,621371 miles)
Step22:  Be aware
Step23:  We dropped about 42% of the data!
Step24:  A little drawing
Step25:  Where is a pickup maximum?
Step26:  The area with the maximum pickups is around the Grand Central Terminal. Not very surprising.
Step27:  So we know that about 97% of the bins have 0 pickups in it. This was expected, because we have static stations only.
Only look at trips in a given bounding box
Step28:  We lost no trips through the bounding box.
Step29:  Let's take a first look at the distribution of the cleaned target variable which we want to estimate
Step30:  Make a new dataframe with features and targets
Step31:  Use minutes for prediction instead of seconds (ceil the time). Definitley more robust than seconds!
Step32:  So we hace 10 different velocities to predict.
Split the data into a training dataset and a test dataset. Evaluate the performance of the decision tree on the test data
Step33:  Changed test size to 5%
Step34:  If you want to export something...
Step35:  Just to be sure
Step36:  So this is a bad result. Let's try random forest now.
Train several random forests and look at their scores
Step37:  Watch at the best results.
Step38:  Train the best random forest again seperately
Step39:  So this plot is not particularly useful in this case. But it's a hinto to how to visualize the results.
Step40:  Dump the best random forest...
<ASSISTANT_TASK:>
Python Code:
data = pd.read_csv('data/bike_oneweekfrom20140505.csv', index_col=0, parse_dates=True)
data.info()
data.columns
Explanation: Use the bash =)
End of explanation
new_column_names = ['trip_time', 'pickup_datetime', 'dropoff_datetime', 'start_station_id',
       'start_station_name', 'pickup_latitude',
       'pickup_longitude', 'end_station_id', 'end_station_name',
       'dropoff_latitude', 'dropoff_longitude', 'bikeid', 'usertype',
       'birth year', 'gender']
data.columns = new_column_names
data.describe()
Explanation: Rename the columns. They're now streamlined with the taxi input.
End of explanation
data['pickup_datetime'] =pd.to_datetime(data['pickup_datetime'], format = '%Y-%m-%d %H:%M:%S')
data['dropoff_datetime'] =pd.to_datetime(data['dropoff_datetime'], format = '%Y-%m-%d %H:%M:%S')
data['trip_time'] = pd.to_timedelta(data['trip_time'], 's')
data.describe().transpose()
data.head()
data.info()
Explanation: Parse the dates.
End of explanation
group_by_start = data.groupby(data.start_station_id, )
usage_freq = group_by_start.trip_time.count()
print('Amount of stations: ' + str(len(usage_freq)))
plt.hist(usage_freq, bins = 50)
plt.title('Frequency of starts from start stations')
Explanation: Have a look at every station. How many pickups are there per station?
End of explanation
group_by_start_name = data.groupby(data.start_station_name, )
usage_freq_name = group_by_start_name.trip_time.count()
usage_freq_name.sort_values().keys()[-10:]
usage_freq_name.sort_values().keys()
Explanation: When you take a look at the map, you won't be surprised that the Grand Central Station is the one with the most frequent usage.
End of explanation
group_by_end_name = data.groupby(data.end_station_name, )
usage_freq_name = group_by_end_name.trip_time.count()
usage_freq_name.sort_values().keys()
Explanation: So we have 326 different start stations.
End of explanation
group_by_start_id = data.groupby(data.start_station_id, )
start_ids = group_by_start_id.trip_time.count().sort_values().keys()
group_by_end_id = data.groupby(data.end_station_id, )
end_ids = group_by_end_id.trip_time.count().sort_values().keys()
len(set(start_ids).intersection(end_ids))
Explanation: We have 326 different end stations, too. But are all end_stations contained in start_stations and vive versa?
End of explanation
data.isnull().sum()
Explanation: So this is the proof, that all stations are used as start and end stations.
Check for missing and false data
End of explanation
(data==0).sum()
Explanation: So there is not that much data missing. That's quite surprising, maybe it's wrong.
End of explanation
from collections import Counter
print(Counter(data.gender))
Counter(data['birth year'])
Explanation: So we only have many zeros in the gender-feature. Gender can have 3 values (0=unknown, 1=male, 2=female)
End of explanation
plt.hist((data['trip_time'] / np.timedelta64(1, 'm')), bins=30, range=[0, 100])
plt.title('Distribution of trip_time in minutes')
plt.xlabel('trip_time')
plt.ylabel('frequency')
plt.savefig('figures/bike_trip_time.png', format='png', dpi=300)
Explanation: So we see, we have also missing values in "gender" and "birth year". But we are not interested in this features, therefore we can ignore this missing values at this point.
Quick pverview about the trip_times
End of explanation
sns.boxplot((data['trip_time'] / np.timedelta64(1, 'm')))
Explanation: We are interested in the trip time in minutes.
End of explanation
data.trip_time.quantile([0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99])
len(data.trip_time.value_counts().values)
Explanation: We have lots of outliers.
So as we can see, we have many outliers.
End of explanation
anomaly = data.loc[(data['dropoff_longitude'].isnull()) | (data['dropoff_latitude'].isnull()) | 
                       (data['pickup_longitude'].isnull()) | (data['pickup_latitude'].isnull())]
data = data.drop(anomaly.index)
anomaly['flag'] = 'geo_NA'
data.isnull().sum()
Explanation: Identify the the cases without geo data and remove them from our data to be processed.
End of explanation
len(data)/(len(data)+len(anomaly))
anomaly.tail()
Explanation: So how many percent of data are left to be processed?
End of explanation
plt.hist(data.trip_time.values / np.timedelta64(1, 'm'), bins=50, range=[0,100])
Explanation: <font color = 'blue'> So we dropped nothing because of missing geo tags. </font color>
End of explanation
data['trip_dist'] = -1 # Init trip_dist
Explanation: We sometimes have some unreasonably small trip_times.
Calculate the distance of a trip. Use the difference of lon / lat data and parse it to metrics.
End of explanation
# inpout for vincenty:(location.latitude, location.longitude)
from geopy.distance import vincenty
for i in range(0, (len(data)-1)):
    pickup = (data.iloc[i]['pickup_latitude'], data.iloc[i]['pickup_longitude'])
    dropoff = (data.iloc[i]['dropoff_latitude'], data.iloc[i]['dropoff_longitude'])
    data.set_value(i, 'trip_dist', vincenty(pickup,dropoff).meters)
data.trip_dist
Explanation: Be aware, this operation is quite slow.
End of explanation
data.columns
Explanation: We cannot use miles because vincenty ceils / floors the results
End of explanation
# data.to_csv('data/bike_20140505_with_dist.csv')
data['avg_velocity'] = data.trip_dist.values/(data.trip_time / (np.timedelta64(1, 'h')))
data.avg_velocity
Explanation: Check if new column is present
End of explanation
data.avg_velocity = data.avg_velocity/1000
data.avg_velocity
Explanation: Parse to km/h
End of explanation
data['avg_velocity'] = data.avg_velocity*0.627371
# data.to_csv('data/bike_20140505_with_dist_and_avg_velo.csv')
plt.hist(data.avg_velocity, bins=60, range=[0,15])
plt.title('Distribution of avg_velocity in mph')
plt.xlabel('avg_velocity')
plt.ylabel('frequency')
plt.savefig('figures/bike_avg_vel.png', format='png', dpi=300)
data.avg_velocity.describe() # in mph
data.head()
data.avg_velocity.quantile([.0001,.01, .5, .75, .95, .975, .99, .995])
Explanation: Convert avg_velocity from km/h to miles/h. (1km = 0,621371 miles)
End of explanation
lb = 5
ub = 15
anomaly = data.loc[(data['avg_velocity'] < lb) | (data['avg_velocity'] > ub)] 
# be careful! Maybe adjust to append.
#anomaly.loc[
anomaly.loc[data.loc[(data['avg_velocity'] < lb)].index,'flag'] = 'too_slow'
anomaly.loc[data.loc[(data['avg_velocity'] > ub)].index,'flag'] = 'too_fast'
data = data.drop(anomaly.index, errors='ignore') # ignore uncontained labels / indices
print(1-len(data)/(len(data)+len(anomaly)))
anomaly.head()
Explanation: Be aware: We calculated the avg_velocity based on the vincenty distance (equals direct connection). A real driver would never drive that fast. We take all avg_velos >= 5mph (approx. >5%-quantile) into account. We use an upper bound, too. It might not be possible to drive "too fast". If someone drives too fast, the technology has had a malfunction.
End of explanation
data.avg_velocity.describe()
anomaly.tail()
Explanation: We dropped about 42% of the data!
End of explanation
x = data.pickup_longitude
y = data.pickup_latitude
bins = 100;
H, xedges, yedges = np.histogram2d(x, y, bins=bins)
plt.jet()
fig = plt.figure(figsize=(20, 10))
#ax.set_title('pcolormesh: exact bin edges')
#mesh = plt.pcolormesh(X, Y, H)
plt.hist2d(x,y,bins=bins)
plt.colorbar()
plt.scatter(xedges[79], yedges[61], marker='x')
ax = fig.gca()
ax.grid(False)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Pickup Frequency')
#ax.set_aspect('equal')
#plt.savefig('figure.pdf', format='pdf')
#plt.savefig('figure.png', format='png')
plt.savefig('figures/bike_stations.png', format='png', dpi=300)
Explanation: A little drawing
End of explanation
np.amax(H)
np.where(H==2001)
Explanation: Where is a pickup maximum?
End of explanation
yedges[79],xedges[61], 
type(H)
print(H.shape)
print(H.size)
print(H.max())
print('Current bin width:')
print(vincenty((xedges[0], yedges[0]), (xedges[1], yedges[1])).meters)
distances = (-1)*np.ones(len(xedges)-1)
for x in range(0,len(xedges)-1, 1):
    distances[x] = vincenty((xedges[x], yedges[x]), (xedges[x+1], yedges[x+1])).meters
print('Sizes of all bins in meters:')
distances
(H==0).sum()/H.size
Explanation: The area with the maximum pickups is around the Grand Central Terminal. Not very surprising.
End of explanation
jfk_geodata = (40.641547, -73.778118)
ridgefield_geodata = (40.856406, -74.020642)
data_in_box = data.loc[(data['dropoff_latitude'] > jfk_geodata[0]) & 
                       (data['dropoff_longitude'] < jfk_geodata[1]) &
                       (data['dropoff_latitude'] < ridgefield_geodata[0]) & 
                       (data['dropoff_longitude'] > ridgefield_geodata[1]) & 
                       (data['pickup_latitude'] > jfk_geodata[0]) & 
                       (data['pickup_longitude'] < jfk_geodata[1]) &
                       (data['pickup_latitude'] < ridgefield_geodata[0]) & 
                       (data['pickup_longitude'] > ridgefield_geodata[1])         
                       ]
# taxidata = taxidata.drop(anomaly.index)
len(data_in_box)/len(data)
Explanation: So we know that about 97% of the bins have 0 pickups in it. This was expected, because we have static stations only.
Only look at trips in a given bounding box
End of explanation
data_in_box.head()
Explanation: We lost no trips through the bounding box.
End of explanation
h = data_in_box.avg_velocity.values
plt.figure(figsize=(20,10))
plt.hist(h, normed=False,  bins=150)
    #, histtype='stepfilled')
#plt.yscale('log')
#plt.ylabel('log(freq x)', fontsize=40)
#plt.xlabel('x = avg_amount_per_minute', fontsize=40)
#print('Min:' +  str(min(h)) + '\nMax:' +  str(max(h)))
#plt.locator_params(axis = 'x', nbins = 20)
plt.show()
data_in_box.head()
Explanation: Let's take a first look at the distribution of the cleaned target variable which we want to estimate:
End of explanation
time_regression_df = pd.DataFrame([ #data_in_box['pickup_datetime'].dt.day, # leave this out
                          data_in_box['pickup_datetime'].dt.dayofweek,
                          data_in_box['pickup_datetime'].dt.hour,
                          data_in_box['pickup_latitude'],
                          data_in_box['pickup_longitude'],
                          data_in_box['dropoff_latitude'],
                          data_in_box['dropoff_longitude'],
                          np.ceil(data_in_box['avg_velocity'])
                         ]).T
time_regression_df.columns = [#'pickup_datetime_day', 
                              'pickup_datetime_dayofweek', 'pickup_datetime_hour',
                                 'pickup_latitude', 'pickup_longitude', 'dropoff_latitude', 'dropoff_longitude',
                                 'avg_velocity_mph']
Explanation: Make a new dataframe with features and targets
End of explanation
time_regression_df.head()
time_regression_df.tail()
time_regression_df.ix[:,0:7].describe()
print(time_regression_df.avg_velocity_mph.value_counts())
print(len(time_regression_df.avg_velocity_mph.value_counts()))
Explanation: Use minutes for prediction instead of seconds (ceil the time). Definitley more robust than seconds!
End of explanation
from sklearn import cross_validation as cv
y = time_regression_df['avg_velocity_mph']
X = time_regression_df.ix[:,0:6]
X_train, X_test, y_train, y_test = cv.train_test_split(X, y,test_size=0.05,random_state=0) 
from sklearn import cross_validation as cv
time_regression_df_train, time_regression_df_test = cv.train_test_split(time_regression_df, test_size=0.05, random_state=99)
y_train = time_regression_df_train['avg_velocity_mph']
x_train = time_regression_df_train.ix[:, 0:6]
y_test = time_regression_df_test['avg_velocity_mph']
x_test = time_regression_df_test.ix[:, 0:6]
Explanation: So we hace 10 different velocities to predict.
Split the data into a training dataset and a test dataset. Evaluate the performance of the decision tree on the test data
End of explanation
x_test.head()
y_test.head()
xy_test = pd.concat([x_test, y_test], axis=1)
xy_test.head()
Explanation: Changed test size to 5%
End of explanation
#Xy_test.to_csv('taxi_tree_test_Xy_20130506-12.csv')
#X_test.to_csv('taxi_tree_test_X_20130506-12.csv')
#y_test.to_csv('taxi_tree_test_y_20130506-12.csv')
# Xy_test.to_csv('bike_tree_test_Xy_20140505-11.csv')
# X_test.to_csv('bike_tree_test_X_20140505-11.csv')
# y_test.to_csv('bike_tree_test_y_20140506-11.csv')
# Xy_test_sample = Xy_test.sample(10000, random_state=99)
# Xy_test_sample.to_csv('taxi_tree_test_Xy_sample.csv')
# Xy_test_sample.head()
Explanation: If you want to export something...
End of explanation
print(x_train.shape)
print(x_train.size)
print(x_test.shape)
print(X.shape)
print(x_train.shape[0]+x_test.shape[0])
import time
# Import the necessary modules and libraries
from sklearn.tree import DecisionTreeRegressor
import numpy as np
import matplotlib.pyplot as plt
regtree = DecisionTreeRegressor(min_samples_split=100, random_state=99, max_depth=20)# formerly 15. 15 is reasonable
                                                                                     # random states: 99
regtree.fit(x_train, y_train)
regtree.score(x_test, y_test)
y_pred = regtree.predict(x_test)
diff_regtree = (y_pred-y_test)
# plt.figure(figsize=(12,10)) # not needed. set values globally
plt.hist(diff_regtree.values, bins=100, range=[-6,6])
print('Perzentile(%): ', [1,5,10,15,25,50,75,90,95,99], '\n', np.percentile(diff_regtree.values, [1,5,10,15,25,50,75,90,95,99]))
print('Absolute time deviation (in 1k): ', sum(abs(diff_regtree))/1000)
plt.title('Simple Decision Tree Regressor')
plt.xlabel('deviation in mph')
plt.ylabel('frequency')
plt.savefig('figures/bike_tree.png', format='png', dpi=300)
diff_regtree.describe()
Explanation: Just to be sure
End of explanation
from sklearn.ensemble import RandomForestRegressor
n_est_list = list(range(2,32,8))
min_sam_leaf_list = [2,50,100,150]
max_depth_list = list(range(2,32,8))
results = np.empty([0,4])
for n_est in n_est_list:
    for min_sam_leaf in min_sam_leaf_list:
        for max_depth in max_depth_list:
            rd_regtree = RandomForestRegressor(n_estimators=n_est,n_jobs=6,min_samples_leaf=min_sam_leaf, random_state=99, max_depth=max_depth)
            rd_regtree.fit(x_train, y_train)
            score = rd_regtree.score(x_test, y_test)
            results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))
best = np.where(results == max(results[:,3]))[0]
results[best,:]
Explanation: So this is a bad result. Let's try random forest now.
Train several random forests and look at their scores
End of explanation
results[np.argsort(results[:, 3])][-10:,:]
results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))
#results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))
#results = np.vstack((results, [n_est, min_sam_leaf, max_depth,score]))
results[:,3]
Explanation: Watch at the best results.
End of explanation
from sklearn.ensemble import RandomForestRegressor
#rd_regtree = RandomForestRegressor(n_estimators=10,n_jobs=6,min_samples_leaf=4, random_state=99, max_depth=20)
rd_regtree = RandomForestRegressor(n_estimators=26,n_jobs=6,min_samples_leaf=2, random_state=99, max_depth=26)
#total sum of diff: 1132
#rd_regtree = RandomForestRegressor(n_estimators=40,n_jobs=-1,min_samples_split=3, random_state=99, max_depth=11)
#total sum of diff: 1129
rd_regtree.fit(x_train, y_train)
print('R²: ', rd_regtree.score(x_test, y_test))
rd_regtree.feature_importances_
tree_list = rd_regtree.estimators_
for i in range(0,len(tree_list)):
    print((tree_list[i].feature_importances_))
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(results[:,0],results[:,2],results[:,3], marker='o')
plt.show()
Explanation: Train the best random forest again seperately
End of explanation
y_pred = rd_regtree.predict(x_test)
np.linalg.norm(np.ceil(y_pred)-y_test)
diff_rd = (y_pred-y_test)
# plt.figure(figsize=(12,10)) # not needed. set values globally
plt.hist(diff_rd.values, bins=100, range=[-6,6])
print('Perzentile(%): ', [1,5,10,15,25,50,75,90,95,99], '\n', np.percentile(diff_rd.values, [1,5,10,15,25,50,75,90,95,99]))
print('Absolute time deviation (in 1k): ', sum(abs(diff_rd))/1000)
plt.title('Random Forest Regressor')
plt.xlabel('deviation in mph')
plt.ylabel('frequency')
plt.savefig('figures/bike_randomforest.png', format='png', dpi=150)
diff_rd.describe()
Explanation: So this plot is not particularly useful in this case. But it's a hinto to how to visualize the results.
End of explanation
from sklearn.externals import joblib
joblib.dump(rd_regtree, 'randforlib/bike_regtree_26x_depth_26_mss_2_PY27.pkl', protocol=2)
! cd randforlib && zip bike_regtree_26x_depth_26_mss_2_PY27.zip bike_regtree_26x_depth_26_mss_2_PY27.pkl*
Explanation: Dump the best random forest...
End of explanation
<END_TASK>