<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Predicting the occupancies of Belgian trains
In this lab, we will go over some of the typical steps in a data science pipeline
Step1:  0. Create a kaggle account! https
Step2:  Processing the data
Cleaning the variables, and adding station variables to our dataframe
Step3:  1.2
Step4:  2.2
Step5:  2.3
Step6:  2.4
Step7:  3. Predictive modeling
Step8:  We train our model on a 'training set' and evaluate it on the testset. Functionality for making this split automatically can be found <a href="http
Step9:  Since we have a lot of 'Null' (+-1/3th) values for our 'class' feature, and we don't want to throw that away, we can try to predict these labels based on the other features, we get +75% accuracy so that seems sufficient. But we can't forgot to do the same thing for the test set!
Step10:  4. 'Advanced' predictive modeling
Step11:  5. Data augmentation with external data sources
There is a unlimited amount of factors that influence the occupancy of a train! Definitely more than the limited amount of data given in the feedback logs. Therefore, we will try to create new features for our dataset using external data sources. Examples of data sources include
Step12:  Transform all null classes to one null class, maybe try to predict the class? Based on to and from and time
Step13:  6. Generating a Kaggle submission and comparing your methodology to others
 6.1
<ASSISTANT_TASK:>
Python Code:
import os
os.getcwd() 
%matplotlib inline
%pylab inline
import pandas as pd
import numpy as np
from collections import Counter, OrderedDict
import json
import matplotlib
import matplotlib.pyplot as plt
import re
from scipy.misc import imread
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import confusion_matrix, accuracy_score, classification_report
Explanation: Predicting the occupancies of Belgian trains
In this lab, we will go over some of the typical steps in a data science pipeline:
Data processing & cleaning
Exploratory Data Analysis
Feature extraction/engineering
Model selection & hyper-parameter tuning
Data linking
...
We will make use of the following technologies and libraries:
Python3.5
Python libraries: pandas, numpy, sklearn, matplotlib, ...
Kaggle
NO SPARK!!! (next lab will deal with machine learning with Spark MLlib)
End of explanation
from pandas.io.json import json_normalize
import pickle
training_json = pd.DataFrame()
with open('data/training_data.nldjson') as data_file:  
        for line in data_file:
            training_json = training_json.append(json_normalize(json.loads(line)))
with open('data/test.nldjson') as data_file:  
        for line in data_file:
            out_test_json = json_normalize(json.loads(line))
out_test = out_test_json
training = training_json
out_test[0:1]
Explanation: 0. Create a kaggle account! https://www.kaggle.com/
The competition can be found here: https://inclass.kaggle.com/c/train-occupancy-prediction-v2/leaderboard
Create an account and form a team (shuffle II), use your names and BDS_ as a prefix in your team name
Note: you can only make 5 submissions per day
There are also student groups from Kortrijk (Master of Science in Industrial Engineering) participating. They get no help at all (you get this notebook) but this is their final lab + they have no project. THEREFORE: Let's push them down the leaderboard!!! ;)  
Your deadline: the end of the kaggle competition. 
Evaluation: Your work will be evaluated for 50%, your result will also matter for another 50%. The top 5 student groups get bonus points for this part of the course!
1. Loading and processing the data
Trains can get really crowded sometimes, so wouldn't it be great to know in advance how busy your train will be, so you can take an earlier or later one? iRail, created just that. their application, SpitsGids, shows you the occupancy of every train in Belgium. Furthermore, you can indicate the occupancy yourself. Using the collected data, machine learning models can be trained to predict what the occupancy level of a train will be.
The dataset which we will use during this labo is composed of two files:
train.nldjson: contains labeled training data (JSON records, separated by newlines)
test.nldjson: unlabeled data for which we will create a submission for a Kaggle competition at the end of this lab (again: JSON records, separated by newlines). Each of the records is uniquely identifiable through an id
A json record has the following structure:
{
 "querytype": "occupancy",
 "querytime": "2016-09-29T16:24:43+02:00",
 "post": {
     "connection": "http://irail.be/connections/008811601/20160929/S85666",
     "from": "http://irail.be/stations/NMBS/008811601",
     "to": "http://irail.be/stations/NMBS/008811676",
     "date": "20160929",
     "vehicle": "http://irail.be/vehicle/S85666",
     "occupancy": "http://api.irail.be/terms/medium"
 },
 "user_agent": "Railer/1610 CFNetwork/808.0.2 Darwin/16.0.0"
}
This is how the five first rows of a processed DataFrame COULD look like
 1.1:  Load in both files and store the data in a pandas DataFrame, different methodologies can be applied in order to parse the JSON records (pd.io.json.json_normalize, json library, ...)
Loading the data
Loading the json files and dumping them via pickle
End of explanation
training['querytime'] = pd.to_datetime(training['querytime'])
out_test['querytime'] = pd.to_datetime(out_test['querytime'])
training = training.dropna()
training['post.occupancy'] = training['post.occupancy'].apply(lambda x: x.split("http://api.irail.be/terms/",1)[1])
training['post.vehicle'] = training['post.vehicle'].apply(lambda x: x.split("http://irail.be/vehicle/",1)[1])
out_test['post.vehicle'] = out_test['post.vehicle'].apply(lambda x: x.split("http://irail.be/vehicle/",1)[1])
#create class column, eg IC058 -> IC
training['post.class'] = training['post.vehicle'].apply(lambda x: " ".join(re.findall("[a-zA-Z]+", x)))
out_test['post.class'] = out_test['post.vehicle'].apply(lambda x: " ".join(re.findall("[a-zA-Z]+", x)))
#reset the index because you have duplicate indexes now because you appended DFs in a for loop
training = training.reset_index()
stations_df = pd.read_csv('data/stations.csv')
stations_df['from'] = stations_df.index
stations_df['destination'] = stations_df['from']
stations_df[0:4]
#post.from en post.to are in the some format of URI
stations_df["zoekterm"]=stations_df["name"]+" trein"
stations_df.loc[stations_df['zoekterm'].str.startswith("Zaventem"), "zoekterm"] = "Zaventem trein"
stations_df.loc[stations_df['zoekterm'].str.startswith("Charleroi"), "zoekterm"] = "Charleroi trein"
stations_df.loc[stations_df['zoekterm'].str.startswith("Brussel"), "zoekterm"] = "Brussel trein"
stations_df.loc[stations_df['zoekterm'].str.startswith("Gent"), "zoekterm"] = "Gent trein"
stations_df.loc[stations_df['zoekterm'].str.startswith("Liège"), "zoekterm"] = "Luik trein"
stations_df.loc[stations_df['zoekterm'].str.startswith("Antwerpen"), "zoekterm"] = "Antwerpen trein" 
druktes_df = pd.read_csv('data/station_druktes.csv')
druktes_df[0:4]
training = pd.merge(training,stations_df[["URI","from"]], left_on = 'post.from', right_on = 'URI')
training = pd.merge(training,stations_df[["URI","destination"]], left_on = 'post.to', right_on = 'URI')
training = training.drop(['URI_y','URI_x'],1)
out_test = pd.merge(out_test,stations_df[["URI","from"]], left_on = 'post.from', right_on = 'URI')
out_test = pd.merge(out_test,stations_df[["URI","destination"]], left_on = 'post.to', right_on = 'URI')
out_test = out_test.drop(['URI_y','URI_x'],1)
Explanation: Processing the data
Cleaning the variables, and adding station variables to our dataframe
End of explanation
fig, ax = plt.subplots(1,1, figsize=(5,5))
training['post.class'].value_counts().plot(kind='pie', ax=ax, autopct='%1.1f%%')
#we have a lot of null/undefined, especially in our test set, we can't simply throw them away
Explanation: 1.2:  Clean the data! Make sure the station- and vehicle-identifiers are in the right format. A station identifier consists of 9 characters (prefix = '00') and a vehicle identifier consists of the concatentation of the vehicle type (IC/L/S/P/...) and the line identifier. Try to fix as much of the records as possible, drop only the unfixable ones. How many records did you drop?
2. Exploratory Data Analysis (EDA)
Let's create some visualisations of our data in order to gain some insights. Which features are useful, which ones aren't?
We will create 3 visualisations:
* Pie chart of the class distribution
* Stacked Bar Chart depicting the distribution for one aggregated variable (such as the weekday or the vehicle type)
* Scattter plot depicting the 'crowdiness' of the stations in Belgium
For each of the visualisations, code to generate the plot has already been handed to you. You only need to prepare the data (i.e. create a new dataframe or select certain columns) such that it complies with the input specifications. If you want to create your own plotting code or extend the given code, you are free to do so!
 2.1:  *Create a pie_chart with the distribution of the different classes. Have a look at our webscraping lab for plotting pie charts. TIP: the value_counts() does most of the work for you!
End of explanation
training['weekday'] = training['querytime'].apply(lambda l: l.weekday()) 
out_test['weekday'] = out_test['querytime'].apply(lambda l: l.weekday()) 
print("timerange from training data:",training['querytime'].min(),training['querytime'].max())
print(training['querytime'].describe())
print(out_test['querytime'].describe())
date_training = training.set_index('querytime')
date_test = out_test.set_index('querytime')
grouper = pd.TimeGrouper("1d")
date_training = date_training.groupby(grouper).size()
date_test = date_test.groupby(grouper).size()
# plot
fig, ax = plt.subplots(1,1, figsize=(10,7))
ax.plot(date_training)
ax.plot(date_test)
fig, ax = plt.subplots(1,1, figsize=(6,6))
training['weekday'].value_counts().plot(kind='pie', ax=ax, autopct='%1.1f%%')
training['post.occupancy'].value_counts()
Explanation: 2.2:  *Analyze the timestamps in the training and testset. First convert the timestamps to a pandas datetime object using pd.datetime. http://pandas.pydata.org/pandas-docs/stable/timeseries.html
Have the column in this data format simplifies a lot of work, since it allows you to convert and extract time features more easily. For example: 
- df['weekday] = df['time'].apply(lambda l: l.weekday())
would map every date to a day of the week in [0,6]. 
A. What are the ranges of training and testset, is your challenges one of interpolating or extrapolating in the future? 
TIP: The describe() function can already be helpful!
B. Plot the number of records in both training and testset per day. Have a look here on how to work with the timegrouper functionality: http://stackoverflow.com/questions/15297053/how-can-i-divide-single-values-of-a-dataframe-by-monthly-averages
C. OPTIONAL: Have insight into the time dependence can get you a long way: Make additional visualizations to make you understand how time affects train occupancy.
End of explanation
training[0:1]
occup = pd.crosstab(training['post.class'], training['post.occupancy'])
weekday = pd.crosstab(training['post.class'], training['weekday'])
occup = occup.drop(['BUS', 'EUR', 'EXT', 'ICE', 'ICT', 'P', 'TGV', 'THA', 'TRN', 'ic', 'null', ''])
occup = occup.apply(lambda r: r/r.sum(), axis=1)
occup[0:4]
weekday = weekday.drop(['BUS', 'EUR', 'EXT', 'ICE', 'ICT', 'P', 'TGV', 'THA', 'TRN', 'ic', 'null', ''])
weekday = weekday.apply(lambda r: r/r.sum(), axis=1)
df_occup = pd.DataFrame(occup)
df_occup.plot.bar(stacked=True);
df_weekday = pd.DataFrame(weekday)
df_weekday.plot.bar(stacked=True);
Explanation: 2.3:  *Create a stacked_bar_chart with the distribution of the three classes over an aggregated variable (group the data by weekday, vehicle_type, ...). More info on creating stacked bar charts can be found here: http://pandas.pydata.org/pandas-docs/stable/visualization.html#bar-plots
The dataframe you need will require your grouping variables as the index, and 1 column occupancy category, for example:
| Index | Occupancy_Low  | Occupancy_Medium | Occupancy_High | Sum_Occupancy |
|-------|----------------|------------------|----------------|---------------|
| IC    | 15             | 30               | 10             | 55 |
| S     | 20             | 10               | 30             | 60 |
| L     | 12             | 9                | 14             | 35 |
If you want the values to be relative (%), add a sum column and use it to divide the occupancy columns
End of explanation
stops = stations_df[['URI','longitude','latitude']] 
dest_count = training['post.to'].value_counts()
dest_count_df = pd.DataFrame({'id':dest_count.index, 'count':dest_count.values})
dest_loc = pd.merge(dest_count_df, stops, left_on = 'id', right_on = 'URI')
dest_loc = dest_loc[['id', 'count', 'latitude','longitude']]
fig, ax = plt.subplots(figsize=(12,10))
ax.scatter(dest_loc.longitude, dest_loc.latitude, s=dest_loc['count'] )
Explanation: 2.4:  * To have an idea about the hotspots in the railway network make a scatter plot that depicts the number of visitors per station. Aggregate on the destination station and use the GTFS dataset at iRail to find the geolocation of the stations (stops.txt): https://gtfs.irail.be/nmbs
End of explanation
def get_seconds_since_midnight(x):
    midnight = x.replace(hour=0, minute=0, second=0, microsecond=0)
    return (x - midnight).seconds
def get_line_number(x):
    pattern = re.compile("^[A-Z]+([0-9]+)$")
    if pattern.match(x):
        return int(pattern.match(x).group(1))
    else:
        return x
    
training['seconds_since_midnight'] = training['querytime'].apply(get_seconds_since_midnight)
training['month'] = training['querytime'].apply(lambda x: x.month)
training['occupancy'] = training['post.occupancy'].map({'low': 0, 'medium': 1, 'high': 2})
out_test['seconds_since_midnight'] = out_test['querytime'].apply(get_seconds_since_midnight)
out_test['month'] = out_test['querytime'].apply(lambda x: x.month)
fig, ax = plt.subplots(figsize=(5, 5))
corr_frame = training[['seconds_since_midnight', 'month', 'occupancy']].corr()
cax = ax.matshow(abs(corr_frame))
fig.colorbar(cax)
tickpos = np.array(range(0,len(corr_frame.columns)))
plt.xticks(tickpos,corr_frame.columns, rotation='vertical')
plt.yticks(tickpos,corr_frame.columns, rotation='horizontal')
plt.grid(None)
pd.tools.plotting.scatter_matrix(training[['seconds_since_midnight', 'month', 'occupancy']],  
                                 alpha=0.2, diagonal='kde', figsize=(10,10))
plt.grid(None)
Explanation: 3. Predictive modeling: creating a baseline
Now that we have processed, cleaned and explored our data it is time to create a predictive model that predicts the occupancies of future Belgian trains. We will start with applying Logistic Regression on features extracted from our initial dataset. Some code has already been given to get you started.
Feature extraction
Some possible features include (bold ones are already implemented for you):
The day of the week
The number of seconds since midnight of the querytime
The train vehicle type (IC/P/L/...)
The line number
The line category
Information about the from- and to-station (their identifier, their coordinates, the number of visitors, ...)
The month
A binary variable indicating whether a morning (6-10AM) or evening jam (3-7PM) is ongoing
...
In order to do reveal relations between these features you can try and plot them with:
<a href="https://datascience.stackexchange.com/questions/10459/calculation-and-visualization-of-correlation-matrix-with-pandas"> Correlation plot </a>
<a href="http://pandas.pydata.org/pandas-docs/stable/visualization.html#visualization-scatter-matrix"> Scatter matrix </a>
These relations can be important since some models do not perform very will when features are highly correlated
Feature normalization
Most models require the features to have a similar range, preferables [0,1]. A minmax scaler is usually sufficient: x -> (x - xmin) / (xmax - xmin)
Scikit will be used quite extensively from now on, have a look here for preprocessing functionality: http://scikit-learn.org/stable/modules/classes.html#module-sklearn.preprocessing
Dealing with categorical variables
All machine learning techniques, except for tree-based methods, assume that variables are ordinal (you can define an order). For some variables, such as the day of the week or the train vehicle type, this is not true. Therefore, a pre-processing step is required that transforms these categorical variables. A few examples of such transformations are:
One-hot-encoding (supported by pandas: get_dummies )
Binary encoding: map each variable to a number, binary encode these numbers and use each bit as a feature (advantage of this technique is that it introduces a lot less new variables in contrast to one-hot-encoding)
Hash encoding
...
 3.1:  Extract more features than the two given ones. Make sure you extract at least one categorical variable, and transform it! What gains (in terms of current accuracy (0.417339475755)) do you achieve with new features in comparison to the given code?
End of explanation
skf = StratifiedKFold(n_splits=5, random_state=1337)
X = training[['seconds_since_midnight', 'month']]
y = training['occupancy']
cms = []
accs = []
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index, :], X.iloc[test_index, :]
    y_train, y_test = y[train_index], y[test_index]
    
    log_reg = LogisticRegression()
    log_reg.fit(X_train, y_train)
    
    predictions = log_reg.predict(X_test)
    cm = confusion_matrix(y_test, predictions)
    cms.append(cm)
    accs.append(accuracy_score(y_test, predictions))
    print(classification_report(y_test, predictions))
    #accs.append(sum([float(cm[i][i]) for i in range(len(cm))])/np.sum(cm))
    
print('Confusion matrix:\n', np.mean(cms, axis=0))
print('Avg accuracy', np.mean(accs), '+-', np.std(accs))
print('Predict all lows', float(len(y[y == 0]))/float(len(y)))
Explanation: We train our model on a 'training set' and evaluate it on the testset. Functionality for making this split automatically can be found <a href="http://scikit-learn.org/stable/modules/classes.html#module-sklearn.model_selection"> here </a>
Our first model is a linear logistic regression model, more information on the API <a href="http://scikit-learn.org/stable/modules/classes.html#module-sklearn.linear_model"> here </a>
The confusion matrix is part of the <a href="http://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics"> metrics functionality </a>
End of explanation
training_class = training_holidays[training_holidays.class_enc != 0]
training_class = training_class[training_class.class_enc != 14]
test_class = training_holidays[(training_holidays.class_enc == 0)|(training_holidays.class_enc == 14)]
training_class["class_pred"]=training_class["class_enc"]
training_holidays_enc = pd.concat([training_class,test_class])
X_train = training_class[['seconds_since_midnight','weekday', 'month','id','id_2']]
X_test = test_class[['seconds_since_midnight','weekday', 'month','id','id_2']]
y_train = training_class['class_enc']
train.occupancy.value_counts()/train.shape[0]
test.occupancy.value_counts()/test.shape[0]
out_test_holidays_druktes.occupancy.value_counts()/out_test_holidays_druktes.shape[0]
from sklearn.linear_model import SGDClassifier
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier
from sklearn.metrics import f1_score
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import VotingClassifier
from sklearn.model_selection import RandomizedSearchCV, GridSearchCV 
from sklearn.cross_validation import train_test_split
train, test = train_test_split(training_holidays_druktes, test_size=0.2, random_state=42)
X_train = train[['seconds_since_midnight','drukte_from','drukte_to','school','name_enc','day__0','day__1','day__2','day__3','day__4','day__5','day__6','from_lat','from_lng','des_lat','des_lng','real_trend','class__0','class__1','class__2','class__3','class__4','class__5','class__6','class__7','class__8', 'temperature', 'humidity', 'windspeed', 'weather_type', 'visibility']]
X_test = test[['seconds_since_midnight','drukte_from','drukte_to','school','name_enc','day__0','day__1','day__2','day__3','day__4','day__5','day__6','from_lat','from_lng','des_lat','des_lng','real_trend','class__0','class__1','class__2','class__3','class__4','class__5','class__6','class__7','class__8', 'temperature', 'humidity', 'windspeed', 'weather_type', 'visibility']]
y_train = train['occupancy']
y_test = test['occupancy']
#month uit de set halen als we ongeziene willen predicten
from xgboost import XGBClassifier
xgb = XGBClassifier(n_estimators=5000, max_depth=3, min_child_weight=6, learning_rate=0.01,
                        colsample_bytree=0.5, subsample=0.6, gamma=0., nthread=-1,
                        max_delta_step=1, objective='multi:softmax')
xgb.fit(X_train, y_train, sample_weight=[1]*len(y_train))
print(xgb.score(X_train,y_train))
print(xgb.score(X_test, y_test))
ac = AdaBoostClassifier()
ada_param_grid = {'n_estimators': [10, 30, 100, 300, 1000],
                  'learning_rate': [0.1, 0.3, 1.0, 3.0]}
ac_grid = GridSearchCV(ac,ada_param_grid,cv=3,
                       scoring='accuracy')
ac_grid.fit(X_train, y_train)
ac = ac_grid.best_estimator_
#ac.fit(X_train, y_train)
#print(ac_grid.score(X_train,y_train))
#print(ac_grid.score(X_test, y_test))
rf = RandomForestClassifier()
param_dist = {"n_estimators": [20],
                "max_depth": [7, None],
              "max_features": range(4, 6),
              "min_samples_split": range(2, 7),
              "min_samples_leaf": range(1, 7),
              "bootstrap": [True, False],
              "criterion": ["gini", "entropy"]}
rand = GridSearchCV(rf,param_dist,cv=3,
                       scoring='accuracy')
rand.fit(X_train, y_train)
rf = rand.best_estimator_
print(rand.best_estimator_)
# rf.fit(X_train, y_train)
# print(rf.score(X_train,y_train))
# print(rf.score(X_test, y_test))
from sklearn.tree import DecisionTreeClassifier
dtc = DecisionTreeClassifier(random_state=0)
dtc.fit(X_train, y_train)
print(dtc.score(X_train,y_train))
print(dtc.score(X_test, y_test))
rf2 = rand.best_estimator_
rf3 = rand.best_estimator_
rf4 = rand.best_estimator_
# voting_clf = VotingClassifier(
# estimators=[('ac', ac), ('rf', rf), ('dtc', dtc),('rf2', rf2), ('rf3', rf3), ('rf4', rf4), ('xgb', xgb)],
# voting='hard'
# )
voting_clf = VotingClassifier(
estimators=[('ac', ac), ('rf', rf), ('xgb', xgb)],
voting='hard'
)
from sklearn.metrics import accuracy_score
for clf in (ac, rf, xgb, voting_clf):
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test)
    print(clf.__class__.__name__, accuracy_score(y_test, y_pred))
ac.fit(X_train, y_train)
voting_clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
print(clf.__class__.__name__, accuracy_score(y_test, y_pred))
pd.DataFrame([X_train.columns, rf.feature_importances_])
y_predict_test = voting_clf.predict(out_test_holidays_druktes[['seconds_since_midnight','drukte_from','drukte_to','school','name_enc','day__0','day__1','day__2','day__3','day__4','day__5','day__6','from_lat','from_lng','des_lat','des_lng','real_trend','class__0','class__1','class__2','class__3','class__4','class__5','class__6','class__7','class__8','temperature', 'humidity', 'windspeed', 'weather_type', 'visibility']])
y_predict_test = rf.predict(out_test_holidays_druktes[['seconds_since_midnight','drukte_from','drukte_to','school','name_enc','day__0','day__1','day__2','day__3','day__4','day__5','day__6','from_lat','from_lng','des_lat','des_lng','real_trend','class__0','class__1','class__2','class__3','class__4','class__5','class__6','class__7','class__8','temperature', 'humidity', 'windspeed', 'weather_type', 'visibility']])
out_test_holidays_druktes["occupancy"] = y_predict_test
out_test_holidays_druktes.occupancy.value_counts()/out_test_holidays_druktes.shape[0]
train.occupancy.value_counts()/train.shape[0]
out_test_holidays_druktes[['seconds_since_midnight','drukte_from','drukte_to','name_enc','class_enc','day__0','day__1','day__2','day__3','day__4','day__5','day__6','trend','occupancy']][0:100]
out_test_holidays_druktes[["id","occupancy"]].to_csv('predictions.csv',index=False)
Explanation: Since we have a lot of 'Null' (+-1/3th) values for our 'class' feature, and we don't want to throw that away, we can try to predict these labels based on the other features, we get +75% accuracy so that seems sufficient. But we can't forgot to do the same thing for the test set!
End of explanation
skf = StratifiedKFold(n_splits=5, random_state=1337)
X = training[['seconds_since_midnight', 'month']]
y = training['occupancy']
cms = []
accs = []
parameters = {#'penalty': ['l1', 'l2'],  # No penalty tuning, cause 'l1' is only supported by liblinear
              # It can be interesting to manually take a look at 'l1' with 'liblinear', since LASSO
              # provides sparse solutions (boils down to the fact that LASSO does some feature selection for you)
              'solver': ['newton-cg', 'lbfgs', 'liblinear', 'sag'],
              'tol': [1e-4, 1e-6, 1e-8], 
              'C': [1e-2, 1e-1, 1.0, 1e1],
              'max_iter': [1e2, 1e3]
             }
for train_index, test_index in skf.split(X, y):
    X_train, X_test = X.iloc[train_index, :], X.iloc[test_index, :]
    y_train, y_test = y[train_index], y[test_index]
    
    tuned_log_reg = GridSearchCV(LogisticRegression(penalty='l2'), parameters, cv=3,
                       scoring='accuracy')
    tuned_log_reg.fit(X_train, y_train)
    print(tuned_log_reg.best_params_)
    
    predictions = tuned_log_reg.predict(X_test)
    cm = confusion_matrix(y_test, predictions)
    cms.append(cm)
    accs.append(accuracy_score(y_test, predictions))
    print(classification_report(y_test, predictions))
    
print('Confusion matrix:\n', np.mean(cms, axis=0))
print('Avg accuracy', np.mean(accs), '+-', np.std(accs))
print('Predict all lows', float(len(y[y == 0]))/float(len(y)))
Explanation: 4. 'Advanced' predictive modeling: model selection & hyper-parameter tuning
Model evaluation and hyper-parameter tuning
In order to evaluate your model, K-fold cross-validation (https://en.wikipedia.org/wiki/Cross-validation_(statistics) ) is often applied. Here, the data is divided in K chunks, K-1 chunks are used for training while 1 chunk is used for testing. Different metrics exist, such as accuracy, AUC, F1 score, and more. For this lab, we will use accuracy.
Some machine learning techniques, supported by sklearn:
SVMs
Decision Trees
Decision Tree Ensemble: AdaBoost, Random Forest, Gradient Boosting
Multi-Level Perceptrons/Neural Networks
Naive Bayes
K-Nearest Neighbor
...
To tune the different hyper-parameters of a machine learning model, again different techniques exist:
* Grid search: exhaustively try all possible parameter combinations (Code to tune the different parameters of our LogReg model has been given)
* Random search: try a number of random combinations, it has been shown that this is quite equivalent to grid search
 4.1:  *Choose one or more machine learning techniques, different from Logistic Regression and apply them to our data, with tuned hyper-parameters! You will see that switching techniques in sklearn is really simple! Which model performs best on this data? *
End of explanation
holiday_pops = pd.read_json('data/holidays.json')
holidays = pd.read_json( (holiday_pops['holidays']).to_json(), orient='index')
holidays['date'] = pd.to_datetime(holidays['date'])
holidays.head(1)
training["date"] = training["querytime"].values.astype('datetime64[D]')
out_test["date"] = out_test["querytime"].values.astype('datetime64[D]')
training_holidays = pd.merge(training,holidays, how="left", on='date')
training_holidays.school = training_holidays.school.fillna(0)
training_holidays.name = training_holidays.name.fillna("geen")
training_holidays[0:1]
out_test_holidays = pd.merge(out_test,holidays, how="left", on='date')
out_test_holidays.school = out_test_holidays.school.fillna(0)
out_test_holidays.name = out_test_holidays.name.fillna("geen")
out_test_holidays[0:1]
from sklearn.preprocessing import LabelEncoder
encoder = LabelEncoder()
#encode the names from the holidays (Summer,Christmas...)
training_holidays["name_enc"] = encoder.fit_transform(training_holidays["name"])
out_test_holidays["name_enc"] = encoder.fit_transform(out_test_holidays["name"])
#encode the classes (IC,TGV,L...)
training_holidays["class_enc"] = encoder.fit_transform(training_holidays["post.class"])
out_test_holidays["class_enc"] = encoder.fit_transform(out_test_holidays["post.class"])
training_holidays=training_holidays.rename(columns = {'too':'destination'})
out_test_holidays=out_test_holidays.rename(columns = {'too':'destination'})
stations_df = pd.merge(stations_df,druktes_df.drop(['Unnamed: 0','station'],1), left_on = 'name', right_on = 'station_link')
Explanation: 5. Data augmentation with external data sources
There is a unlimited amount of factors that influence the occupancy of a train! Definitely more than the limited amount of data given in the feedback logs. Therefore, we will try to create new features for our dataset using external data sources. Examples of data sources include:
Weather APIs
A holiday calendar
Event calendars
Connection and delay information of the SpitsGidsAPI
Data from the NMBS/SNCB
Twitter and other social media
many, many more
In order to save time, a few 'prepared' files have already been given to you. Of course, you are free to scrape/generate your own data as well:
Hourly weather data for all stations in Belgium, from August till April weather_data.zip
A file which contains the vehicle identifiers and the stations where this vehicle stops  line_info.csv
Based on this line_info, you can construct a graph of the rail net in Belgium and apply some fancy graph features (pagerank, edge betweenness, ...) iGraph experiments.ipynb
A file containing the coordinates of a station, and the number of visitors during week/weekend for 2015  station_druktes.csv
A file with some of the holidays (this can definitely be extended) holidays.json
For event data, there is the Eventful API
 5.1:  Pick one (or more) external data source(s) and link your current data frame to that data source (requires some creativity in most cases). Extract features from your new, linked data source and re-train your model. How much gain did you achieve?
Als we kijken naar "training.id.value_counts()" dan zien we vooral dat het om studenten bestemmingen gaat, misschien komt dat omdat het vooral hen zijn die deze app gebruiken? We moeten dus nadenken wanneer zij de trein nemen, en wat dat kan beinvloeden. Misschien het aantal studenten per station incorporeren?
End of explanation
def transform_druktes(row):
    start = row['from']
    destination = row['destination']
    day = row['weekday']
    row['from_lat']=stations_df[stations_df["from"] == start]["latitude"].values[0]
    row['from_lng']=stations_df[stations_df["destination"] == destination]["longitude"].values[0]
    row['des_lat']=stations_df[stations_df["from"] == start]["latitude"].values[0]
    row['des_lng']=stations_df[stations_df["destination"] == destination]["longitude"].values[0]
    row['zoekterm']=stations_df[stations_df["destination"] == destination]["zoekterm"].values[0]
    if day == 5:
        row['drukte_from']=stations_df[stations_df["from"] == start]["zaterdag"].values[0]
        row['drukte_to']=stations_df[stations_df["destination"] == destination]["zaterdag"].values[0]
    elif day == 6:
        row['drukte_from']=stations_df[stations_df["from"] == start]["zondag"].values[0]
        row['drukte_to']=stations_df[stations_df["destination"] == destination]["zondag"].values[0]
    elif day == 4:
        row['drukte_from']=stations_df[stations_df["from"] == start]["zondag"].values[0]/5.0*1.11
        row['drukte_to']=stations_df[stations_df["destination"] == destination]["zondag"].values[0]/5.0*1.11
    elif day == 3:
        row['drukte_from']=stations_df[stations_df["from"] == start]["zondag"].values[0]/5.0*1.21
        row['drukte_to']=stations_df[stations_df["destination"] == destination]["zondag"].values[0]/5.0*1.21
    elif day == 2:
        row['drukte_from']=stations_df[stations_df["from"] == start]["zondag"].values[0]/5.0*0.736
        row['drukte_to']=stations_df[stations_df["destination"] == destination]["zondag"].values[0]/5.0*0.736
    elif day == 1:
        row['drukte_from']=stations_df[stations_df["from"] == start]["zondag"].values[0]/5.0*0.92
        row['drukte_to']=stations_df[stations_df["destination"] == destination]["zondag"].values[0]/5.*0.92
    else:
        row['drukte_from']=stations_df[stations_df["from"] == start]["week"].values[0]/5.0*1.016
        row['drukte_to']=stations_df[stations_df["destination"] == destination]["week"].values[0]/5.0*1.016
    
    return row
training_holidays_druktes = training_holidays_druktes.apply(transform_druktes, axis=1)
out_test_holidays_druktes = training_holidays_druktes.apply(transform_druktes, axis=1)
training_holidays_druktes = pd.concat([training_holidays_druktes,
           pd.get_dummies(training_holidays_druktes['weekday'], prefix="day_"), 
           ],1)
out_test_holidays_druktes = pd.concat([out_test_holidays_druktes,
           pd.get_dummies(out_test_holidays_druktes['weekday'], prefix="day_"), 
           ],1)
trends_df = pd.DataFrame()
real_trend_df = pd.DataFrame()
from pytrends.request import TrendReq
import pandas as pd
# enter your own credentials
google_username = "davidjohansmolders@gmail.com"
google_password = "*******"
#path = ""
# Login to Google. Only need to run this once, the rest of requests will use the same session.
pytrend = TrendReq(google_username, google_password, custom_useragent='My Pytrends Script')
for i in range(0,645):
    if i % 4 != 0:
        continue
    try:
        pytrend.build_payload(kw_list=[stations_df[stations_df.destination == i].zoekterm.values[0], stations_df[stations_df.destination == i+1].zoekterm.values[0], stations_df[stations_df.destination == i+2].zoekterm.values[0], stations_df[stations_df.destination == i+3].zoekterm.values[0], "Brussel trein"],geo="BE",timeframe='2016-07-27 2017-04-05')
        real_trend_df = pd.concat([real_trend_df,pytrend.interest_over_time()], axis=1)
    except:
        continue
no_dup_trends = trends_df.T.groupby(level=0).first().T
training_holidays_druktes = pd.merge(training_holidays_druktes,stations_df[["destination","zoekterm"]], left_on = 'destination', right_on = 'destination')
out_test_holidays_druktes = pd.merge(out_test_holidays_druktes,stations_df[["destination","zoekterm"]], left_on = 'destination', right_on = 'destination')
int(real_trend_df.loc["2016-07-28"]["Brussel trein"])
training_holidays_druktes_copy = training_holidays_druktes
out_test_holidays_druktes_copy = out_test_holidays_druktes
training_holidays_druktes = training_holidays_druktes_copy
out_test_holidays_druktes = out_test_holidays_druktes_copy
def get_trends(row):
    zoek = str(row.zoekterm) 
    datum = str(row["date"])[0:10]
    try:
        row["real_trend"] = int(real_trend_df.loc[datum][zoek])
    except:
        row["real_trend"] = 0
    return row
training_holidays_druktes = training_holidays_druktes.apply(get_trends, axis=1)
out_test_holidays_druktes = out_test_holidays_druktes.apply(get_trends, axis=1)
training_holidays_druktes = training_holidays_druktes.drop(['post.date','post.from','post.vehicle','querytype','user_agent','post.to','name','post.class'],1)
out_test_holidays_druktes = out_test_holidays_druktes.drop(['post.date','post.from','post.vehicle','querytype','user_agent','post.to','name','post.class'],1)
training_holidays_druktes["hour"] = training_holidays_druktes["querytime"].values.astype('datetime64[h]').astype('str')
out_test_holidays_druktes["hour"] = out_test_holidays_druktes["querytime"].values.astype('datetime64[h]').astype('str')
training_holidays_druktes["hour_lag"] = (training_holidays_druktes["querytime"].values.astype('datetime64[h]')-2).astype('str')
out_test_holidays_druktes["hour_lag"] = (out_test_holidays_druktes["querytime"].values.astype('datetime64[h]')-2).astype('str')
training_holidays_druktes["timeframe"] = training_holidays_druktes["hour_lag"]+" "+training_holidays_druktes["hour"]
out_test_holidays_druktes["timeframe"] = out_test_holidays_druktes["hour_lag"]+" "+out_test_holidays_druktes["hour"]
# enter your own credentials
google_username = "davidjohansmolders@gmail.com"
google_password = "*******"
#path = ""
# Login to Google. Only need to run this once, the rest of requests will use the same session.
pytrend = TrendReq(google_username, google_password, custom_useragent='My Pytrends Script')
def get_hour_trends(row):
    zoek = str(row.zoekterm_x) 
    tijd = str(row["timeframe"])
    try:
        pytrend.build_payload(kw_list=[zoek],timeframe=tijd)
        row["hour_trend"] = int(pytrend.interest_over_time()[zoek].sum())
    except:
        row["hour_trend"] = 0
    return row
training_holidays_druktes = training_holidays_druktes.apply(get_hour_trends, axis=1)
out_test_holidays_druktes = out_test_holidays_druktes.apply(get_hour_trends, axis=1)
training_holidays_druktes = pd.concat([training_holidays_druktes,
           pd.get_dummies(training_holidays_druktes['class_enc'], prefix="class_"), 
           ],1)
out_test_holidays_druktes = pd.concat([out_test_holidays_druktes,
           pd.get_dummies(out_test_holidays_druktes['class_enc'], prefix="class_"), 
           ],1)
# file names for all csv files containing weather information per month
weather_csv = ['weather_data_apr_1', 'weather_data_apr_2', 'weather_data_aug_1', 'weather_data_aug_2', 'weather_data_dec_1', 'weather_data_dec_2', 'weather_data_feb_1', 'weather_data_feb_2', 'weather_data_jan_1', 'weather_data_jan_2', 'weather_data_july_1', 'weather_data_july_2', 'weather_data_mar_1', 'weather_data_mar_2', 'weather_data_nov_1', 'weather_data_nov_2', 'weather_data_oct_1', 'weather_data_oct_2', 'weather_data_sep_1', 'weather_data_sep_2']
for i in range(len(weather_csv)):
    weather_csv[i] = 'data/weather_data/' + weather_csv[i] + '.csv'
# create column of station index
stations_df['station_index'] = stations_df.index
# put all weather data in an array
weather_months = []
for csv in weather_csv:
    weather_month = pd.read_csv(csv)
    
    # convert date_time to a datetime object
    weather_month['date_time'] = pd.to_datetime(weather_month['date_time'])
    weather_month = weather_month.drop(['Unnamed: 0','lat','lng'], 1)
    
    weather_months.append(weather_month)
    
# concatenate all weather data
weather = pd.concat(weather_months)
# merge weather month with station to convert station name to index (that can be found in holiday_druktes)
weather = pd.merge(weather, stations_df[["name", "station_index"]], left_on = 'station_name', right_on = 'name')
weather = weather.drop(['station_name', 'name'], 1)
# truncate querytime to the hour in new column
training_holidays_druktes['querytime_hour'] = training_holidays_druktes['querytime'].apply(lambda dt: datetime.datetime(dt.year, dt.month, dt.day, dt.hour))
# join training with weather data
training_holidays_druktes_weather = pd.merge(training_holidays_druktes, weather, how='left', left_on = ['destination', 'querytime_hour'], right_on = ['station_index', 'date_time'])
training_holidays_druktes_weather = training_holidays_druktes_weather.drop(['querytime_hour', 'date_time', 'station_index'], 1)
training_holidays_druktes_weather = training_holidays_druktes_weather.drop_duplicates()
# fill null rows of weather data with there mean
training_holidays_druktes_weather['temperature'].fillna(training_holidays_druktes_weather['temperature'].mean(), inplace=True)
training_holidays_druktes_weather['humidity'].fillna(training_holidays_druktes_weather['humidity'].mean(), inplace=True)
training_holidays_druktes_weather['windspeed'].fillna(training_holidays_druktes_weather['windspeed'].mean(), inplace=True)
training_holidays_druktes_weather['visibility'].fillna(training_holidays_druktes_weather['visibility'].mean(), inplace=True)
training_holidays_druktes_weather['weather_type'].fillna(training_holidays_druktes_weather['weather_type'].mean(), inplace=True)
#cast weather type to int
training_holidays_druktes_weather['weather_type'] = training_holidays_druktes_weather['weather_type'].astype(int)
# Add weather data to test data
# truncate querytime to the hour in new column
out_test_holidays_druktes['querytime_hour'] = out_test_holidays_druktes['querytime'].apply(lambda dt: datetime.datetime(dt.year, dt.month, dt.day, dt.hour))
# join test with weather data
out_test_holidays_druktes_weather = pd.merge(out_test_holidays_druktes, weather, how='left', left_on = ['destination', 'querytime_hour'], right_on = ['station_index', 'date_time'])
out_test_holidays_druktes_weather = out_test_holidays_druktes_weather.drop(['querytime_hour', 'date_time', 'station_index'], 1)
out_test_holidays_druktes_weather = out_test_holidays_druktes_weather.drop_duplicates()
# fill null rows of weather data with there mean
out_test_holidays_druktes_weather['temperature'].fillna(out_test_holidays_druktes_weather['temperature'].mean(), inplace=True)
out_test_holidays_druktes_weather['humidity'].fillna(out_test_holidays_druktes_weather['humidity'].mean(), inplace=True)
out_test_holidays_druktes_weather['windspeed'].fillna(out_test_holidays_druktes_weather['windspeed'].mean(), inplace=True)
out_test_holidays_druktes_weather['visibility'].fillna(out_test_holidays_druktes_weather['visibility'].mean(), inplace=True)
out_test_holidays_druktes_weather['weather_type'].fillna(out_test_holidays_druktes_weather['weather_type'].mean(), inplace=True)
#cast weather type to int
out_test_holidays_druktes_weather['weather_type'] = out_test_holidays_druktes_weather['weather_type'].astype(int)
# set out_test_holidays_druktes and training_holidays_druktes equal to weather counterpart such that we don't need to change all variable names above
out_test_holidays_druktes = out_test_holidays_druktes_weather
training_holidays_druktes = training_holidays_druktes_weather
Explanation: Transform all null classes to one null class, maybe try to predict the class? Based on to and from and time
End of explanation
pickle.dump(training_holidays_druktes,open("temp_data/training_holidays_druktes.pkl","wb"))
pickle.dump(out_test_holidays_druktes,open("temp_data/out_test_holidays_druktes.pkl","wb"))
training_holidays_druktes = pd.read_pickle("temp_data/training_holidays_druktes.pkl")
out_test_holidays_druktes = pd.read_pickle("temp_data/out_test_holidays_druktes.pkl")
training_holidays_druktes[0:5]
Explanation: 6. Generating a Kaggle submission and comparing your methodology to others
 6.1:  Train your best performing model on train.nldjson and generate predictions for test.nldjson. Create a file called submission.csv with format listed below and submit it on the Kaggle competition! What's your position on the leaderboard?
End of explanation
<END_TASK>