Kaggle San Francisco Crime Classification

Berkeley MIDS W207 Final Project: Sam Goodgame, Sarah Cha, Kalvin Kao, Bryan Moore

Basic Modeling

Environment and Data


In [1]:
# Import relevant libraries:
import time
import numpy as np
import pandas as pd
from sklearn.neighbors import KNeighborsClassifier
from sklearn import preprocessing
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import StandardScaler
from sklearn.naive_bayes import BernoulliNB
from sklearn.naive_bayes import MultinomialNB
from sklearn.naive_bayes import GaussianNB
from sklearn.grid_search import GridSearchCV
from sklearn.metrics import classification_report
from sklearn.metrics import log_loss
from sklearn.linear_model import LogisticRegression
from sklearn import svm
from sklearn.neural_network import MLPClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
# Import Meta-estimators
from sklearn.ensemble import AdaBoostClassifier
from sklearn.ensemble import BaggingClassifier
from sklearn.ensemble import GradientBoostingClassifier
# Import Calibration tools
from sklearn.calibration import CalibratedClassifierCV

# Set random seed and format print output:
np.random.seed(0)
np.set_printoptions(precision=3)


C:\Users\kalvi\AppData\Local\conda\conda\envs\py2All1\lib\site-packages\sklearn\cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.
  "This module will be removed in 0.20.", DeprecationWarning)
C:\Users\kalvi\AppData\Local\conda\conda\envs\py2All1\lib\site-packages\sklearn\grid_search.py:43: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. This module will be removed in 0.20.
  DeprecationWarning)

DDL to construct table for SQL transformations:

CREATE TABLE kaggle_sf_crime (
dates TIMESTAMP,                                
category VARCHAR,
descript VARCHAR,
dayofweek VARCHAR,
pd_district VARCHAR,
resolution VARCHAR,
addr VARCHAR,
X FLOAT,
Y FLOAT);

Getting training data into a locally hosted PostgreSQL database:

\copy kaggle_sf_crime FROM '/Users/Goodgame/Desktop/MIDS/207/final/sf_crime_train.csv' DELIMITER ',' CSV HEADER;

SQL Query used for transformations:

SELECT
  category,
  date_part('hour', dates) AS hour_of_day,
  CASE
    WHEN dayofweek = 'Monday' then 1
    WHEN dayofweek = 'Tuesday' THEN 2
    WHEN dayofweek = 'Wednesday' THEN 3
    WHEN dayofweek = 'Thursday' THEN 4
    WHEN dayofweek = 'Friday' THEN 5
    WHEN dayofweek = 'Saturday' THEN 6
    WHEN dayofweek = 'Sunday' THEN 7
  END AS dayofweek_numeric,
  X,
  Y,
  CASE
    WHEN pd_district = 'BAYVIEW' THEN 1
    ELSE 0
  END AS bayview_binary,
    CASE
    WHEN pd_district = 'INGLESIDE' THEN 1
    ELSE 0
  END AS ingleside_binary,
    CASE
    WHEN pd_district = 'NORTHERN' THEN 1
    ELSE 0
  END AS northern_binary,
    CASE
    WHEN pd_district = 'CENTRAL' THEN 1
    ELSE 0
  END AS central_binary,
    CASE
    WHEN pd_district = 'BAYVIEW' THEN 1
    ELSE 0
  END AS pd_bayview_binary,
    CASE
    WHEN pd_district = 'MISSION' THEN 1
    ELSE 0
  END AS mission_binary,
    CASE
    WHEN pd_district = 'SOUTHERN' THEN 1
    ELSE 0
  END AS southern_binary,
    CASE
    WHEN pd_district = 'TENDERLOIN' THEN 1
    ELSE 0
  END AS tenderloin_binary,
    CASE
    WHEN pd_district = 'PARK' THEN 1
    ELSE 0
  END AS park_binary,
    CASE
    WHEN pd_district = 'RICHMOND' THEN 1
    ELSE 0
  END AS richmond_binary,
    CASE
    WHEN pd_district = 'TARAVAL' THEN 1
    ELSE 0
  END AS taraval_binary
FROM kaggle_sf_crime;

Load the data into training, development, and test:


In [2]:
data_path = "./data/train_transformed.csv"

df = pd.read_csv(data_path, header=0)
x_data = df.drop('category', 1)
y = df.category.as_matrix()

# Impute missing values with mean values:
x_complete = x_data.fillna(x_data.mean())
X_raw = x_complete.as_matrix()

# Scale the data between 0 and 1:
X = MinMaxScaler().fit_transform(X_raw)

# Shuffle data to remove any underlying pattern that may exist:
shuffle = np.random.permutation(np.arange(X.shape[0]))
X, y = X[shuffle], y[shuffle]

# Separate training, dev, and test data:
test_data, test_labels = X[800000:], y[800000:]
dev_data, dev_labels = X[700000:800000], y[700000:800000]
train_data, train_labels = X[:700000], y[:700000]

mini_train_data, mini_train_labels = X[:75000], y[:75000]
mini_dev_data, mini_dev_labels = X[75000:100000], y[75000:100000]
labels_set = set(mini_dev_labels)
print(labels_set)
print(len(labels_set))


set(['KIDNAPPING', 'WEAPON LAWS', 'SECONDARY CODES', 'WARRANTS', 'PROSTITUTION', 'EMBEZZLEMENT', 'LOITERING', 'SUICIDE', 'DRIVING UNDER THE INFLUENCE', 'VEHICLE THEFT', 'ROBBERY', 'BURGLARY', 'SUSPICIOUS OCC', 'ARSON', 'BRIBERY', 'FORGERY/COUNTERFEITING', 'BAD CHECKS', 'DRUNKENNESS', 'GAMBLING', 'OTHER OFFENSES', 'FRAUD', 'RECOVERED VEHICLE', 'DRUG/NARCOTIC', 'TRESPASS', 'LARCENY/THEFT', 'VANDALISM', 'MISSING PERSON', 'EXTORTION', 'LIQUOR LAWS', 'SEX OFFENSES NON FORCIBLE', 'SEX OFFENSES FORCIBLE', 'STOLEN PROPERTY', 'ASSAULT', 'FAMILY OFFENSES', 'NON-CRIMINAL', 'DISORDERLY CONDUCT', 'RUNAWAY'])
37

Loading the data, version 2, with weather features to improve performance:

We seek to add features to our models that will improve performance with respect to out desired performance metric. There is evidence that there is a correlation between weather patterns and crime, with some experts even arguing for a causal relationship between weather and crime [1]. More specifically, a 2013 paper published in Science showed that higher temperatures and extreme rainfall led to large increases in conflict. In the setting of strong evidence that weather influences crime, we see it as a candidate for additional features to improve the performance of our classifiers. Weather data was gathered from (insert source). Certain features from this data set were incorporated into the original crime data set in order to add features that were hypothesizzed to improve performance. These features included (insert what we eventually include).


In [6]:
data_path = "./data/train_transformed.csv"

df = pd.read_csv(data_path, header=0)
x_data = df.drop('category', 1)
y = df.category.as_matrix()

########## Adding the date back into the data
import csv
import time
import calendar
data_path = "./data/train.csv"
dataCSV = open(data_path, 'rt')
csvData = list(csv.reader(dataCSV))
csvFields = csvData[0] #['Dates', 'Category', 'Descript', 'DayOfWeek', 'PdDistrict', 'Resolution', 'Address', 'X', 'Y']
allData = csvData[1:]
dataCSV.close()

df2 = pd.DataFrame(allData)
df2.columns = csvFields
dates = df2['Dates']
dates = dates.apply(time.strptime, args=("%Y-%m-%d %H:%M:%S",))
dates = dates.apply(calendar.timegm)
#print(dates.head())

x_data['secondsFromEpoch'] = dates
colnames = x_data.columns.tolist()
colnames = colnames[-1:] + colnames[:-1]
x_data = x_data[colnames]
##########

########## Adding the weather data into the original crime data
weatherData1 = "./data/1027175.csv"
weatherData2 = "./data/1027176.csv"
dataCSV = open(weatherData1, 'rt')
csvData = list(csv.reader(dataCSV))
csvFields = csvData[0] #['Dates', 'Category', 'Descript', 'DayOfWeek', 'PdDistrict', 'Resolution', 'Address', 'X', 'Y']
allWeatherData1 = csvData[1:]
dataCSV.close()

dataCSV = open(weatherData2, 'rt')
csvData = list(csv.reader(dataCSV))
csvFields = csvData[0] #['Dates', 'Category', 'Descript', 'DayOfWeek', 'PdDistrict', 'Resolution', 'Address', 'X', 'Y']
allWeatherData2 = csvData[1:]
dataCSV.close()

weatherDF1 = pd.DataFrame(allWeatherData1)
weatherDF1.columns = csvFields
dates1 = weatherDF1['DATE']
sunrise1 = weatherDF1['DAILYSunrise']
sunset1 = weatherDF1['DAILYSunset']

weatherDF2 = pd.DataFrame(allWeatherData2)
weatherDF2.columns = csvFields
dates2 = weatherDF2['DATE']
sunrise2 = weatherDF2['DAILYSunrise']
sunset2 = weatherDF2['DAILYSunset']

#functions for processing the sunrise and sunset times of each day
def get_hour_and_minute(milTime):
    hour = int(milTime[:-2])
    minute = int(milTime[-2:])
    return [hour, minute]

def get_date_only(date):
    return time.struct_time(tuple([date[0], date[1], date[2], 0, 0, 0, date[6], date[7], date[8]]))

def structure_sun_time(timeSeries, dateSeries):
    sunTimes = timeSeries.copy()
    for index in range(len(dateSeries)):
        sunTimes[index] = time.struct_time(tuple([dateSeries[index][0], dateSeries[index][1], dateSeries[index][2], timeSeries[index][0], timeSeries[index][1], dateSeries[index][5], dateSeries[index][6], dateSeries[index][7], dateSeries[index][8]]))
    return sunTimes

dates1 = dates1.apply(time.strptime, args=("%Y-%m-%d %H:%M",))
sunrise1 = sunrise1.apply(get_hour_and_minute)
sunrise1 = structure_sun_time(sunrise1, dates1)
sunrise1 = sunrise1.apply(calendar.timegm)
sunset1 = sunset1.apply(get_hour_and_minute)
sunset1 = structure_sun_time(sunset1, dates1)
sunset1 = sunset1.apply(calendar.timegm)
dates1 = dates1.apply(calendar.timegm)

dates2 = dates2.apply(time.strptime, args=("%Y-%m-%d %H:%M",))
sunrise2 = sunrise2.apply(get_hour_and_minute)
sunrise2 = structure_sun_time(sunrise2, dates2)
sunrise2 = sunrise2.apply(calendar.timegm)
sunset2 = sunset2.apply(get_hour_and_minute)
sunset2 = structure_sun_time(sunset2, dates2)
sunset2 = sunset2.apply(calendar.timegm)
dates2 = dates2.apply(calendar.timegm)

weatherDF1['DATE'] = dates1
weatherDF1['DAILYSunrise'] = sunrise1
weatherDF1['DAILYSunset'] = sunset1
weatherDF2['DATE'] = dates2
weatherDF2['DAILYSunrise'] = sunrise2
weatherDF2['DAILYSunset'] = sunset2

weatherDF = pd.concat([weatherDF1,weatherDF2[32:]],ignore_index=True)

# Starting off with some of the easier features to work with-- more to come here . . . still in beta
weatherMetrics = weatherDF[['DATE','HOURLYDRYBULBTEMPF','HOURLYRelativeHumidity', 'HOURLYWindSpeed', \
                            'HOURLYSeaLevelPressure', 'HOURLYVISIBILITY', 'DAILYSunrise', 'DAILYSunset']]
weatherMetrics = weatherMetrics.convert_objects(convert_numeric=True)
weatherDates = weatherMetrics['DATE']
#'DATE','HOURLYDRYBULBTEMPF','HOURLYRelativeHumidity', 'HOURLYWindSpeed',
#'HOURLYSeaLevelPressure', 'HOURLYVISIBILITY'
timeWindow = 10800 #3 hours
hourlyDryBulbTemp = []
hourlyRelativeHumidity = []
hourlyWindSpeed = []
hourlySeaLevelPressure = []
hourlyVisibility = []
dailySunrise = []
dailySunset = []
daylight = []
test = 0
for timePoint in dates:#dates is the epoch time from the kaggle data
    relevantWeather = weatherMetrics[(weatherDates <= timePoint) & (weatherDates > timePoint - timeWindow)]
    hourlyDryBulbTemp.append(relevantWeather['HOURLYDRYBULBTEMPF'].mean())
    hourlyRelativeHumidity.append(relevantWeather['HOURLYRelativeHumidity'].mean())
    hourlyWindSpeed.append(relevantWeather['HOURLYWindSpeed'].mean())
    hourlySeaLevelPressure.append(relevantWeather['HOURLYSeaLevelPressure'].mean())
    hourlyVisibility.append(relevantWeather['HOURLYVISIBILITY'].mean())
    dailySunrise.append(relevantWeather['DAILYSunrise'].iloc[-1])
    dailySunset.append(relevantWeather['DAILYSunset'].iloc[-1])
    daylight.append(1.0*((timePoint >= relevantWeather['DAILYSunrise'].iloc[-1]) and (timePoint < relevantWeather['DAILYSunset'].iloc[-1])))
    #if timePoint < relevantWeather['DAILYSunset'][-1]:
        #daylight.append(1)
    #else:
        #daylight.append(0)
    
    if test%100000 == 0:
        print(relevantWeather)
    test += 1

hourlyDryBulbTemp = pd.Series.from_array(np.array(hourlyDryBulbTemp))
hourlyRelativeHumidity = pd.Series.from_array(np.array(hourlyRelativeHumidity))
hourlyWindSpeed = pd.Series.from_array(np.array(hourlyWindSpeed))
hourlySeaLevelPressure = pd.Series.from_array(np.array(hourlySeaLevelPressure))
hourlyVisibility = pd.Series.from_array(np.array(hourlyVisibility))
dailySunrise = pd.Series.from_array(np.array(dailySunrise))
dailySunset = pd.Series.from_array(np.array(dailySunset))
daylight = pd.Series.from_array(np.array(daylight))

x_data['HOURLYDRYBULBTEMPF'] = hourlyDryBulbTemp
x_data['HOURLYRelativeHumidity'] = hourlyRelativeHumidity
x_data['HOURLYWindSpeed'] = hourlyWindSpeed
x_data['HOURLYSeaLevelPressure'] = hourlySeaLevelPressure
x_data['HOURLYVISIBILITY'] = hourlyVisibility
x_data['DAILYSunrise'] = dailySunrise
x_data['DAILYSunset'] = dailySunset
x_data['Daylight'] = daylight

x_data.to_csv(path_or_buf="C:/MIDS/W207 final project/x_data.csv")
##########

# Impute missing values with mean values:
x_complete = x_data.fillna(x_data.mean())
X_raw = x_complete.as_matrix()

# Scale the data between 0 and 1:
X = MinMaxScaler().fit_transform(X_raw)

# Shuffle data to remove any underlying pattern that may exist:
shuffle = np.random.permutation(np.arange(X.shape[0]))
X, y = X[shuffle], y[shuffle]

# Separate training, dev, and test data:
test_data, test_labels = X[800000:], y[800000:]
dev_data, dev_labels = X[700000:800000], y[700000:800000]
train_data, train_labels = X[:700000], y[:700000]

mini_train_data, mini_train_labels = X[:75000], y[:75000]
mini_dev_data, mini_dev_labels = X[75000:100000], y[75000:100000]
labels_set = set(mini_dev_labels)
print(labels_set)
print(len(labels_set))
print(train_data[:10])


C:\Users\kalvi\AppData\Local\conda\conda\envs\py2All1\lib\site-packages\ipykernel_launcher.py:102: FutureWarning: convert_objects is deprecated.  Use the data-type specific converters pd.to_datetime, pd.to_timedelta and pd.to_numeric.
              DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
140303  1431550560                56.0                    70.0   
140304  1431554040                57.0                    67.0   
140305  1431554160                57.0                    69.0   
140306  1431554400                57.0                    69.0   
140307  1431557760                57.0                    74.0   
140308  1431558900                57.0                    72.0   

        HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
140303              0.0                   29.93              10.0   
140304              0.0                     NaN              10.0   
140305              8.0                   29.92              10.0   
140306              8.0                   29.92               NaN   
140307              9.0                   29.91              10.0   
140308              8.0                     NaN              10.0   

        DAILYSunrise  DAILYSunset  
140303    1431493260   1431544260  
140304    1431493260   1431544260  
140305    1431493260   1431544260  
140306    1431493260   1431544260  
140307    1431493260   1431544260  
140308    1431493260   1431544260  
              DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
124567  1389459360                54.0                    77.0   
124568  1389462960                53.0                    80.0   
124569  1389464640                54.0                    77.0   
124570  1389465540                52.0                    82.0   
124571  1389466560                52.0                    83.0   

        HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
124567             14.0                   30.22              10.0   
124568              9.0                   30.23              10.0   
124569             10.0                     NaN              10.0   
124570             11.0                     NaN              10.0   
124571              7.0                   30.25              10.0   

        DAILYSunrise  DAILYSunset  
124567    1389425040   1389460320  
124568    1389425040   1389460320  
124569    1389425040   1389460320  
124570    1389425040   1389460320  
124571    1389425040   1389460320  
              DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
109394  1348037760                56.0                    77.0   
109395  1348041360                57.0                    74.0   
109396  1348044960                57.0                    74.0   

        HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
109394             14.0                   30.12              10.0   
109395              9.0                   30.13              10.0   
109396             15.0                   30.14              10.0   

        DAILYSunrise  DAILYSunset  
109394    1348034100   1348078200  
109395    1348034100   1348078200  
109396    1348034100   1348078200  
             DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
92913  1301896560                50.0                    66.0   
92914  1301900160                52.0                    64.0   
92915  1301903760                56.0                    55.0   

       HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
92913              0.0                   30.17              10.0   
92914              0.0                   30.17              10.0   
92915              0.0                   30.17              10.0   

       DAILYSunrise  DAILYSunset  
92913    1301896260   1301942100  
92914    1301896260   1301942100  
92915    1301896260   1301942100  
             DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
77034  1254243360                60.0                    60.0   
77035  1254246960                58.0                    67.0   
77036  1254250560                57.0                    72.0   

       HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
77034             39.0                   29.99              10.0   
77035             37.0                   29.99              10.0   
77036             25.0                   30.01              10.0   

       DAILYSunrise  DAILYSunset  
77034    1254204240   1254246840  
77035    1254204240   1254246840  
77036    1254204240   1254246840  
             DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
61113  1209480960                58.0                    60.0   
61114  1209484560                57.0                    58.0   
61115  1209484800                57.0                    58.0   
61116  1209488160                55.0                    64.0   

       HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
61113             32.0                   30.07             10.00   
61114             33.0                   30.06             10.00   
61115             33.0                   30.06              9.94   
61116             34.0                   30.05             10.00   

       DAILYSunrise  DAILYSunset  
61113    1209446100   1209495540  
61114    1209446100   1209495540  
61115    1209446040   1209495540  
61116    1209446040   1209495540  
             DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
43577  1162742160                65.0                    73.0   
43578  1162742400                65.0                    73.0   
43579  1162745760                63.0                    75.0   
43580  1162749360                62.0                    78.0   

       HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
43577             15.0                   30.20             10.00   
43578             15.0                   30.20              9.94   
43579             14.0                   30.19             10.00   
43580             14.0                   30.21             10.00   

       DAILYSunrise  DAILYSunset  
43577    1162708740   1162746420  
43578    1162708800   1162746360  
43579    1162708800   1162746360  
43580    1162708800   1162746360  
             DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
27369  1118696400                56.0                    81.0   
27370  1118699760                55.0                    82.0   
27371  1118700000                55.0                    83.0   
27372  1118703600                55.0                    83.0   

       HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
27369             17.0                   29.86              8.95   
27370             13.0                   29.85             10.00   
27371             13.0                   29.85              9.94   
27372             10.0                   29.85              9.94   

       DAILYSunrise  DAILYSunset  
27369    1118638020   1118691120  
27370    1118638020   1118691120  
27371    1118638020   1118691120  
27372    1118638020   1118691120  
             DATE  HOURLYDRYBULBTEMPF  HOURLYRelativeHumidity  \
11993  1074722400                51.0                    77.0   
11994  1074726000                48.0                    68.0   
11995  1074729600                49.0                    71.0   

       HOURLYWindSpeed  HOURLYSeaLevelPressure  HOURLYVISIBILITY  \
11993              0.0                   30.29               NaN   
11994              3.0                   30.30              9.94   
11995              0.0                   30.31               NaN   

       DAILYSunrise  DAILYSunset  
11993    1074669600   1074705720  
11994    1074669600   1074705720  
11995    1074756000   1074792120  
set(['KIDNAPPING', 'WEAPON LAWS', 'SECONDARY CODES', 'WARRANTS', 'PROSTITUTION', 'EMBEZZLEMENT', 'LOITERING', 'FRAUD', 'DRIVING UNDER THE INFLUENCE', 'SEX OFFENSES FORCIBLE', 'ROBBERY', 'BURGLARY', 'SUSPICIOUS OCC', 'FAMILY OFFENSES', 'BRIBERY', 'FORGERY/COUNTERFEITING', 'BAD CHECKS', 'DRUNKENNESS', 'GAMBLING', 'OTHER OFFENSES', 'RECOVERED VEHICLE', 'SUICIDE', 'ARSON', 'DRUG/NARCOTIC', 'TRESPASS', 'LARCENY/THEFT', 'VANDALISM', 'NON-CRIMINAL', 'EXTORTION', 'PORNOGRAPHY/OBSCENE MAT', 'LIQUOR LAWS', 'SEX OFFENSES NON FORCIBLE', 'VEHICLE THEFT', 'STOLEN PROPERTY', 'ASSAULT', 'MISSING PERSON', 'DISORDERLY CONDUCT', 'RUNAWAY'])
38
[[  7.986e-01   3.913e-01   6.667e-01   5.101e-02   1.450e-03   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
    1.000e+00   0.000e+00   0.000e+00   0.000e+00   3.911e-01   7.627e-01
    7.186e-02   4.993e-01   9.000e-01   7.987e-01   7.987e-01   1.000e+00]
 [  6.125e-01   1.000e+00   6.667e-01   5.611e-02   2.273e-04   0.000e+00
    1.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   3.799e-01   8.192e-01
    2.455e-01   4.913e-01   9.735e-01   6.124e-01   6.124e-01   0.000e+00]
 [  9.715e-01   2.174e-01   0.000e+00   4.993e-02   1.101e-03   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   1.000e+00   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   1.830e-01   8.249e-01
    8.982e-02   8.420e-01   1.000e+00   9.716e-01   9.716e-01   0.000e+00]
 [  5.213e-01   6.522e-01   1.000e+00   4.562e-02   1.440e-03   0.000e+00
    0.000e+00   1.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   5.810e-01   4.162e-01
    3.752e-01   5.108e-01   1.000e+00   5.213e-01   5.213e-01   1.000e+00]
 [  4.211e-01   7.391e-01   3.333e-01   5.120e-02   1.205e-03   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   1.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   3.665e-01   7.220e-01
    3.784e-01   6.753e-01   9.988e-01   4.211e-01   4.211e-01   1.000e+00]
 [  5.076e-01   9.565e-01   0.000e+00   4.213e-02   1.232e-03   0.000e+00
    0.000e+00   1.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   2.584e-01   6.893e-01
    4.790e-01   5.108e-01   9.985e-01   5.075e-01   5.075e-01   0.000e+00]
 [  7.153e-01   5.217e-01   1.000e+00   5.048e-02   1.475e-03   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
    1.000e+00   0.000e+00   0.000e+00   0.000e+00   3.925e-01   5.932e-01
    2.635e-01   5.606e-01   9.985e-01   7.153e-01   7.153e-01   1.000e+00]
 [  3.559e-01   5.217e-01   3.333e-01   4.487e-02   6.184e-04   0.000e+00
    1.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   4.511e-01   5.650e-01
    3.653e-01   5.909e-01   9.985e-01   3.559e-01   3.559e-01   1.000e+00]
 [  7.701e-01   8.261e-01   1.667e-01   1.188e-02   4.986e-04   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   1.000e+00   4.721e-01   6.017e-01
    5.269e-01   3.810e-01   1.000e+00   7.701e-01   7.701e-01   0.000e+00]
 [  4.562e-01   7.826e-01   1.000e+00   5.475e-02   1.292e-03   0.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   0.000e+00   1.000e+00
    0.000e+00   0.000e+00   0.000e+00   0.000e+00   5.852e-01   6.271e-01
    6.168e-01   3.810e-01   9.985e-01   4.561e-01   4.561e-01   1.000e+00]]

Formatting to meet Kaggle submission standards


In [42]:
# The Kaggle submission format requires listing the ID of each example.
# This is to remember the order of the IDs after shuffling
allIDs = np.array(list(df.axes[0]))
allIDs = allIDs[shuffle]

testIDs = allIDs[800000:]
devIDs = allIDs[700000:800000]
trainIDs = allIDs[:700000]

# Extract the column names for the required submission format
sampleSubmission_path = "./data/sampleSubmission.csv"
sampleDF = pd.read_csv(sampleSubmission_path)
allColumns = list(sampleDF.columns)
featureColumns = allColumns[1:]

# Extracting the test data for a baseline submission
real_test_path = "./data/test_transformed.csv"
testDF = pd.read_csv(real_test_path, header=0)
real_test_data = testDF

test_complete = real_test_data.fillna(real_test_data.mean())
Test_raw = test_complete.as_matrix()

TestData = MinMaxScaler().fit_transform(Test_raw)

# Here we remember the ID of each test data point, in case we ever decide to shuffle the test data for some reason
testIDs = list(testDF.axes[0])

Generate baseline prediction probabilities from MNB classifier and store in a .csv file


In [43]:
# Generate a baseline MNB classifier and make it return prediction probabilities for the actual test data
def MNB():
    mnb = MultinomialNB(alpha = 0.0000001)
    mnb.fit(train_data, train_labels)
    #print("\n\nMultinomialNB accuracy on dev data:", mnb.score(dev_data, dev_labels))
    return mnb.predict_proba(real_test_data)
MNB()

baselinePredictionProbabilities = MNB()

# Place the resulting prediction probabilities in a .csv file in the required format
# First, turn the prediction probabilties into a data frame
resultDF = pd.DataFrame(baselinePredictionProbabilities,columns=featureColumns)
# Add the IDs as a final column
resultDF.loc[:,'Id'] = pd.Series(testIDs,index=resultDF.index)
# Make the 'Id' column the first column
colnames = resultDF.columns.tolist()
colnames = colnames[-1:] + colnames[:-1]
resultDF = resultDF[colnames]
# Output to a .csv file
resultDF.to_csv('result.csv',index=False)

Note: the code above will shuffle data differently every time it's run, so model accuracies will vary accordingly.


In [44]:
## Data sub-setting quality check-point
print(train_data[:1])
print(train_labels[:1])


[[  7.391e-01   6.667e-01   6.117e-02   9.047e-04   1.000e+00   0.000e+00
    0.000e+00   0.000e+00   1.000e+00   0.000e+00   0.000e+00   0.000e+00
    0.000e+00   0.000e+00   0.000e+00]]
['WEAPON LAWS']

In [45]:
# Modeling quality check-point with MNB--fast model

def MNB():
    mnb = MultinomialNB(alpha = 0.0000001)
    mnb.fit(train_data, train_labels)
    print("\n\nMultinomialNB accuracy on dev data:", mnb.score(dev_data, dev_labels))
    
MNB()



MultinomialNB accuracy on dev data: 0.22314

Defining Performance Criteria

As determined by the Kaggle submission guidelines, the performance criteria metric for the San Francisco Crime Classification competition is Multi-class Logarithmic Loss (also known as cross-entropy). There are various other performance metrics that are appropriate for different domains: accuracy, F-score, Lift, ROC Area, average precision, precision/recall break-even point, and squared error.

(Describe each performance metric and a domain in which it is preferred. Give Pros/Cons if able)

  • Multi-class Log Loss:

  • Accuracy:

  • F-score:

  • Lift:

  • ROC Area:

  • Average precision

  • Precision/Recall break-even point:

  • Squared-error:

Model Prototyping

We will start our classifier and feature engineering process by looking at the performance of various classifiers with default parameter settings in predicting labels on the mini_dev_data:


In [46]:
def model_prototype(train_data, train_labels, eval_data, eval_labels):
    knn = KNeighborsClassifier(n_neighbors=5).fit(train_data, train_labels)
    bnb = BernoulliNB(alpha=1, binarize = 0.5).fit(train_data, train_labels)
    mnb = MultinomialNB().fit(train_data, train_labels)
    log_reg = LogisticRegression().fit(train_data, train_labels)
    support_vm = svm.SVC().fit(train_data, train_labels)
    neural_net = MLPClassifier().fit(train_data, train_labels)
    random_forest = RandomForestClassifier().fit(train_data, train_labels)
    decision_tree = DecisionTreeClassifier().fit(train_data, train_labels)
    
    models = [knn, bnb, mnb, log_reg, support_vm, neural_net, random_forest, decision_tree]
    for model in models:
        eval_prediction_probabilities = model.predict_proba(eval_data)
        eval_predictions = model.predict(eval_data)
        set_eval_predictions = set(eval_predictions)
        crime_labels = list(set_eval_predictions)
        print(model, "Multi-class Log Loss:", log_loss(y_true = eval_labels, y_pred = eval_prediction_probabilities, labels = crime_labels), "\n\n")

# model_prototype(mini_train_data, mini_train_labels, mini_dev_data, mini_dev_labels)

Adding Features, Hyperparameter Tuning, and Model Calibration To Improve Prediction For Each Classifier

Here we seek to optimize the performance of our classifiers in a three-step, dynamnic engineering process.

1) Feature addition

We previously added components from the weather data into the original SF crime data as new features. Here, we will incorporate those features into our classifiers to determine whether or not they improve performance as hypothesized.

2) Hyperparameter tuning

Each classifier has parameters that we can engineer to further optimize performance, as opposed to using the default parameter values as we did above in the model prototyping cell. (Will likely use pipeline here in future)

3) Model calibration

We can calibrate the models via Platt Scaling or Isotonic Regression to attempt to improve their performance.

  • Platt Scaling: (brief explanation of how it works)

  • Isotonic Regression: ((brief explanation of how it works))

Feature Addition

(Decide whether to do this as one big step, or to engineer the addition of individual features for individual classifiers. Can likely do all feature addition here, continuing on the work done with the weather data import completed in earlier cells)

K-Nearest Neighbors

Feature addition:
Hyperparameter tuning:

For the KNN classifier, we can seek to optimize the following classifier parameters: n-neighbors, weights, and the power parameter ('p').


In [ ]:
# def k_neighbors_tuned():
    
    k_value_tuning = [i for i in range(1,100,2)]
    weight_tuning = ['uniform', 'distance']
    power_parameter_tuning = [1,2]
    for k in k_value_tuning:
        for w in weight_tuning:
            for p in power_parameter_tuning:
                tuned_KNN = KNeighborsClassifier(n_neighbors=k, weights=w, p=p).fit(train_data, train_labels)
                dev_preds = tuned_KNN.predict(dev_data)
                set_dev_preds = set(dev_preds)
                crime_labels_tuned = list(set_dev_preds)
                dev_prediction_probabilities = tuned_KNN.predict_proba(dev_data)
                print(tuned_KNN, "Multi-class Log Loss:", log_loss(y_true = dev_labels, y_pred = dev_prediction_probabilities, labels = crime_labels_tuned), "\n\n")

# k_neighbors_tuned()
Model calibration:

We will consider embeding this step within the for loop for the hyperparameter tuning. More likely we will pipeline it along with the hyperparameter tuning steps. We will then use GridSearchCV top find the optimized parameters based on our performance metric of Mutli-Class Log Loss.


In [ ]:
# def k_calibrated():

# Here we will calibrate the KNN classifier with both Platt Scaling and with Isotonic Regression using CalibratedClassifierCV
# with various parameter settings.  The "method" parameter can be set to "sigmoid" or to "isotonic", 
# corresponding to Platt Scaling and to Isotonic Regression respectively.

    methods = ['signoid', 'isotonic']
    for m in methods:
        CCV_for_KNN = CalibratedClassifierCV()
        
        # Will likely embed this step within the for loop for the hyperparameter tuning as that makes more sense. #
        # Or will pipeline it along with the hyperparameter tuning steps #

# k_calibrated()

Multinomial, Bernoulli, and Gaussian Naive Bayes


In [ ]:
def GNB():
    gnb = GaussianNB()
    gnb.fit(train_data, train_labels)
    print("GaussianNB accuracy on dev data:", 
          gnb.score(dev_data, dev_labels))
    
    # Gaussian Naive Bayes requires the data to have a relative normal distribution. Sometimes
    # adding noise can improve performance by making the data more normal:
    train_data_noise = np.random.rand(train_data.shape[0],train_data.shape[1])
    modified_train_data = np.multiply(train_data,train_data_noise)    
    gnb_noise = GaussianNB()
    gnb.fit(modified_train_data, train_labels)
    print("GaussianNB accuracy with added noise:", 
          gnb.score(dev_data, dev_labels))    
    
# Going slightly deeper with hyperparameter tuning and model calibration:
def BNB(alphas):
    
    bnb_one = BernoulliNB(binarize = 0.5)
    bnb_one.fit(train_data, train_labels)
    print("\n\nBernoulli Naive Bayes accuracy when alpha = 1 (the default value):",
          bnb_one.score(dev_data, dev_labels))
    
    bnb_zero = BernoulliNB(binarize = 0.5, alpha=0)
    bnb_zero.fit(train_data, train_labels)
    print("BNB accuracy when alpha = 0:", bnb_zero.score(dev_data, dev_labels))
    
    bnb = BernoulliNB(binarize=0.5)
    clf = GridSearchCV(bnb, param_grid = alphas)
    clf.fit(train_data, train_labels)
    print("Best parameter for BNB on the dev data:", clf.best_params_)
    
    clf_tuned = BernoulliNB(binarize = 0.5, alpha=0.00000000000000000000001)
    clf_tuned.fit(train_data, train_labels)
    print("Accuracy using the tuned Laplace smoothing parameter:", 
          clf_tuned.score(dev_data, dev_labels), "\n\n")
    

def investigate_model_calibration(buckets, correct, total):
    clf_tuned = BernoulliNB(binarize = 0.5, alpha=0.00000000000000000000001)
    clf_tuned.fit(train_data, train_labels)
    
    # Establish data sets
    pred_probs = clf_tuned.predict_proba(dev_data)
    max_pred_probs = np.array(pred_probs.max(axis=1))
    preds = clf_tuned.predict(dev_data)
        
    # For each bucket, look at the predictions that the model yields. 
    # Keep track of total & correct predictions within each bucket.
    bucket_bottom = 0
    bucket_top = 0
    for bucket_index, bucket in enumerate(buckets):
        bucket_top = bucket
        for pred_index, pred in enumerate(preds):
            if (max_pred_probs[pred_index] <= bucket_top) and (max_pred_probs[pred_index] > bucket_bottom):
                total[bucket_index] += 1
                if preds[pred_index] == dev_labels[pred_index]:
                    correct[bucket_index] += 1
        bucket_bottom = bucket_top

def MNB():
    mnb = MultinomialNB(alpha = 0.0000001)
    mnb.fit(train_data, train_labels)
    print("\n\nMultinomialNB accuracy on dev data:", mnb.score(dev_data, dev_labels))

alphas = {'alpha': [0.00000000000000000000001, 0.0000001, 0.0001, 0.001, 
                    0.01, 0.1, 0.0, 0.5, 1.0, 2.0, 10.0]}
buckets = [0.5, 0.9, 0.99, 0.999, .9999, 0.99999, 1.0]
correct = [0 for i in buckets]
total = [0 for i in buckets]

MNB()
GNB()
BNB(alphas)
investigate_model_calibration(buckets, correct, total)

for i in range(len(buckets)):
   accuracy = 0.0
   if (total[i] > 0): accuracy = correct[i] / total[i]
   print('p(pred) <= %.13f    total = %3d    accuracy = %.3f' %(buckets[i], total[i], accuracy))

The Bernoulli Naive Bayes and Multinomial Naive Bayes models can predict whether a loan will be good or bad with XXX% accuracy.

Hyperparameter tuning:

We will prune the work above. Will seek to optimize the alpha parameter (Laplace smoothing parameter) for MNB and BNB classifiers.

Model calibration:

Here we will calibrate the MNB, BNB and GNB classifiers with both Platt Scaling and with Isotonic Regression using CalibratedClassifierCV with various parameter settings. The "method" parameter can be set to "sigmoid" or to "isotonic", corresponding to Platt Scaling and to Isotonic Regression respectively. Will likely embed this step within the for loop for the hyperparameter tuning as that makes more sense. Or will pipeline it along with the hyperparameter tuning steps. We will then use GridSearchCV top find the optimized parameters based on our performance metric of Mutli-Class Log Loss.

THE REST OF THE MODEL CALIBRATION SECTIONS ARE SIMILAR AND THE OUTLINE WILL NOT BE REPEATED AS WOULD BE REDUNDANT.

Logistic Regression

Hyperparameter tuning:

For the Logistic Regression classifier, we can seek to optimize the following classifier parameters: penalty (l1 or l2), C (inverse of regularization strength), solver ('newton-cg', 'lbfgs', 'liblinear', or 'sag')

Model calibration:

See above

Decision Tree

Hyperparameter tuning:

For the Decision Tree classifier, we can seek to optimize the following classifier parameters: min_samples_leaf (the minimum number of samples required to be at a leaf node), max_depth

From readings, setting min_samples_leaf to approximately 1% of the data points can stop the tree from inappropriately classifying outliers, which can help to improve accuracy (unsure if significantly improves MCLL).

Model calibration:

See above

Support Vector Machines

Hyperparameter tuning:

For the SVM classifier, we can seek to optimize the following classifier parameters: C (penalty parameter C of the error term), kernel ('linear', 'poly', 'rbf', sigmoid', or 'precomputed')

See source [2] for parameter optimization in SVM

Model calibration:

See above

Neural Nets

Hyperparameter tuning:

For the Neural Networks MLP classifier, we can seek to optimize the following classifier parameters: hidden_layer_sizes, activation ('identity', 'logistic', 'tanh', 'relu'), solver ('lbfgs','sgd', adam'), alpha, learning_rate ('constant', 'invscaling','adaptive')

Model calibration:

See above

Random Forest

Hyperparameter tuning:

For the Random Forest classifier, we can seek to optimize the following classifier parameters: n_estimators (the number of trees in the forsest), max_features, max_depth, min_samples_leaf, bootstrap (whether or not bootstrap samples are used when building trees), oob_score (whether or not out-of-bag samples are used to estimate the generalization accuracy)

Model calibration:

See above

Meta-estimators

AdaBoost Classifier

Hyperparameter tuning:

There are no major changes that we seek to make in the AdaBoostClassifier with respect to default parameter values.

Adaboosting each classifier:

We will run the AdaBoostClassifier on each different classifier from above, using the classifier settings with optimized Multi-class Log Loss after hyperparameter tuning and calibration.

Bagging Classifier

Hyperparameter tuning:

For the Bagging meta classifier, we can seek to optimize the following classifier parameters: n_estimators (the number of trees in the forsest), max_samples, max_features, bootstrap (whether or not bootstrap samples are used when building trees), bootstrap_features (whether features are drawn with replacement), and oob_score (whether or not out-of-bag samples are used to estimate the generalization accuracy)

Bagging each classifier:

We will run the BaggingClassifier on each different classifier from above, using the classifier settings with optimized Multi-class Log Loss after hyperparameter tuning and calibration.

Gradient Boosting Classifier

Hyperparameter tuning:

For the Gradient Boosting meta classifier, we can seek to optimize the following classifier parameters: n_estimators (the number of trees in the forsest), max_depth, min_samples_leaf, and max_features

Gradient Boosting each classifier:

We will run the GradientBoostingClassifier with loss = 'deviance' (as loss = 'exponential' uses the AdaBoost algorithm) on each different classifier from above, using the classifier settings with optimized Multi-class Log Loss after hyperparameter tuning and calibration.

Final evaluation on test data


In [ ]:
# Here we will likely use Pipeline and GridSearchCV in order to find the overall classifier with optimized Multi-class Log Loss.
# This will be the last step after all attempts at feature addition, hyperparameter tuning, and calibration are completed
# and the corresponding performance metrics are gathered.

References

1) Hsiang, Solomon M. and Burke, Marshall and Miguel, Edward. "Quantifying the Influence of Climate on Human Conflict". Science, Vol 341, Issue 6151, 2013

2) Huang, Cheng-Lung. Wang, Chieh-Jen. "A GA-based feature selection and parameters optimization for support vector machines". Expert Systems with Applications, Vol 31, 2006, p 231-240

3) More to come


In [ ]: