<jupyter_start><jupyter_text>Assignment - DS Unit 2  / Sprint 4 /  Module 3  Loading and Feature Engineering<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.model_selection import train_test_split

pd.options.display.max_columns = 500
pd.options.display.max_rows = 500
# Download and join data into a dataframe
def load(): 
    fremont_bridge = 'https://data.seattle.gov/api/views/65db-xm6k/rows.csv?accessType=DOWNLOAD'
    
    bicycle_weather = 'https://raw.githubusercontent.com/jakevdp/PythonDataScienceHandbook/master/notebooks/data/BicycleWeather.csv'

    counts = pd.read_csv(fremont_bridge, index_col='Date', parse_dates=True, 
                         infer_datetime_format=True)

    weather = pd.read_csv(bicycle_weather, index_col='DATE', parse_dates=True, 
                          infer_datetime_format=True)

    daily = counts.resample('d').sum()
    daily['Total'] = daily.sum(axis=1)
    daily = daily[['Total']] # remove other columns

    weather_columns = ['PRCP', 'SNOW', 'SNWD', 'TMAX', 'TMIN', 'AWND']
    daily = daily.join(weather[weather_columns], how='inner')
    
    # Make a feature for yesterday's total
    daily['Total_yesterday'] = daily.Total.shift(1)
    daily = daily.drop(index=daily.index[0])
    
    return daily

daily = load()
# Modified from code cells 17-21 at
# https://jakevdp.github.io/PythonDataScienceHandbook/05.06-linear-regression.html#Example:-Predicting-Bicycle-Traffic

def jake_wrangle(X):  
    X = X.copy()

    # patterns of use generally vary from day to day; 
    # let's add binary columns that indicate the day of the week:
    days = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
    for i, day in enumerate(days):
        X[day] = (X.index.dayofweek == i).astype(float)


    # we might expect riders to behave differently on holidays; 
    # let's add an indicator of this as well:
    from pandas.tseries.holiday import USFederalHolidayCalendar
    cal = USFederalHolidayCalendar()
    holidays = cal.holidays('2012', '2016')
    X = X.join(pd.Series(1, index=holidays, name='holiday'))
    X['holiday'].fillna(0, inplace=True)


    # We also might suspect that the hours of daylight would affect 
    # how many people ride; let's use the standard astronomical calculation 
    # to add this information:
    def hours_of_daylight(date, axis=23.44, latitude=47.61):
        """Compute the hours of daylight for the given date"""
        days = (date - pd.datetime(2000, 12, 21)).days
        m = (1. - np.tan(np.radians(latitude))
             * np.tan(np.radians(axis) * np.cos(days * 2 * np.pi / 365.25)))
        return 24. * np.degrees(np.arccos(1 - np.clip(m, 0, 2))) / 180.

    X['daylight_hrs'] = list(map(hours_of_daylight, X.index))

    
    # temperatures are in 1/10 deg C; convert to C
    X['TMIN'] /= 10
    X['TMAX'] /= 10
    
    # We can also calcuate the average temperature.
    X['Temp (C)'] = 0.5 * (X['TMIN'] + X['TMAX'])

    # precip is in 1/10 mm; convert to inches
    X['PRCP'] /= 254

    # In addition to the inches of precipitation, let's add a flag that 
    # indicates whether a day is dry (has zero precipitation):
    X['dry day'] = (X['PRCP'] == 0).astype(int)


    # Let's add a counter that increases from day 1, and measures how many 
    # years have passed. This will let us measure any observed annual increase 
    # or decrease in daily crossings:
    X['annual'] = (X.index - X.index[0]).days / 365.

    return X

#X_train = jake_wrangle(X_train)
#X_test = jake_wrangle(X_test)
daily = jake_wrangle(daily)
# light rain = 0 - 0.10 PRCP
# Moderate rain = 0.10 - 0.30 PRCP
# Heavy rain = 0.30+
daily['light_rain'] = (daily['PRCP'] > 0) & (daily['PRCP'] < 0.10)
daily['moderate_rain'] = (daily['PRCP'] >= 0.1) & (daily['PRCP'] < 0.30)
daily['heavy_rain'] = (daily['PRCP'] >= 0.30)
daily['weekend_day'] = (daily['Sat']==1) | (daily['Sun'] == 1)
train = daily[:-100] # everything but the last 100 days
test  = daily[-100:] # last 100 days
train.shape, test.shape
X_train = train.drop(columns='Total')
y_train = train['Total']

X_test  = test.drop(columns='Total')
y_test  = test['Total']

X_train.shape, y_train.shape, X_test.shape, y_test.shape<jupyter_output><empty_output><jupyter_text>Fitting Data to Model<jupyter_code>from sklearn.model_selection import RandomizedSearchCV
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
from sklearn.metrics import mean_absolute_error
params= {
    'n_estimators': [50, 100, 150, 200],
    'max_depth': [3, 4, 5, 6, None],
    'criterion' : ['mae'],
     
}

search = RandomizedSearchCV(
    RandomForestRegressor(random_state=42),
    param_distributions = params,
    scoring='neg_mean_absolute_error',
    n_iter = 20,
    cv=3,
    n_jobs = -1,
    verbose=True,
    return_train_score=True
)

search.fit(X_train, y_train)<jupyter_output>Fitting 3 folds for each of 20 candidates, totalling 60 fits<jupyter_text>Scoring Feature Importances<jupyter_code>search.best_score_
best = search.best_estimator_
n = len(X_train.columns)
figsize = (5, 10)

importances = pd.Series(best.feature_importances_, X_train.columns)
top_n = importances.sort_values()[-n:]
plt.figure(figsize = figsize)
top_n.plot.barh(color='gray');<jupyter_output><empty_output><jupyter_text>Permutation Importance<jupyter_code>!pip install eli5
import eli5
from eli5.sklearn import PermutationImportance

permuter = PermutationImportance(best, cv='prefit',
                                n_iter=2, random_state=42)

permuter.fit(X_test.values, y_test)
feature_names = X_test.columns.tolist()
eli5.show_weights(permuter, top=None, feature_names=feature_names)<jupyter_output><empty_output><jupyter_text>Partial Dependence Plots<jupyter_code># Install pdpbox
!pip install pdpbox
# Partial Dependence Plot
from pdpbox.pdp import pdp_isolate, pdp_plot

feature = 'Total_yesterday'

isolated = pdp_isolate(
    model = best,
    dataset=X_test,
    model_features=X_test.columns,
    feature = feature

)

pdp_plot(isolated, feature_name=feature);
from pdpbox.pdp import pdp_interact, pdp_interact_plot

features = ['Total_yesterday', 'PRCP']


interaction = pdp_interact(
    model = best,
    dataset=X_test,
    model_features = X_test.columns,
    features = features
)

pdp_interact_plot(interaction, plot_type='grid', feature_names = features, figsize = (15,8));<jupyter_output><empty_output><jupyter_text>Shapley Values<jupyter_code>!pip install shap
# Check for date to predict
X_test.head()
data_for_prediction = X_test[X_test.index=='2015-05-30']
data_for_prediction
# Shapley Plot
import shap
shap.initjs()
explainer = shap.TreeExplainer(best)
shap_values = explainer.shap_values(data_for_prediction)
shap.force_plot(explainer.expected_value, shap_values, data_for_prediction)
X_test.tail()
data_for_prediction = X_test[X_test.index=='2015-05-28']
data_for_prediction
shap.initjs()
explainer = shap.TreeExplainer(best)
shap_values = explainer.shap_values(data_for_prediction)
shap.force_plot(explainer.expected_value, shap_values, data_for_prediction)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1247009148486981, 'fraction_numerical': 0.025193525686136524, 'mean_word_length': 4.267605633802817, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 27, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '153817', 'license': 'MIT', 'n_tokens_mistral': 2703, 'n_tokens_neox': 2469, 'n_words': 567, 'path': 'Connor_Heraty_Model_Interpretation_Assignment.ipynb', 'repo_name': 'connorpheraty/DS-Unit-2-Sprint-4-Practicing-Understanding', 'size': '7105'}