# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ![texto alternativo](https://drive.google.com/uc?id=1CXVh-Dy6T792nWexE0kH45AOqqFj2_wG)
#
#
# <center> <h1>Maestría en Analytics</h1> </center>
# <center> <h1>Modelos Avanzados para el Análisis de Datos 1</h1> </center>
# <center> <h1>Exercise 7 - Car Price Prediction </h1> </center>
#
# **Grupo 3**
# - <NAME>
# - <NAME> 
# - <NAME>
# - <NAME>
#
# 29 June 2020

# # Exercise 7
#
# # Part 1 - DT
#
# ## Capital Bikeshare data

# ## Introduction
#
# - Capital Bikeshare dataset from Kaggle: [data](https://github.com/justmarkham/DAT8/blob/master/data/bikeshare.csv), [data dictionary](https://www.kaggle.com/c/bike-sharing-demand/data)
# - Each observation represents the bikeshare rentals initiated during a given hour of a given day

# %matplotlib inline
import pandas as pd
import numpy as np
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor, export_graphviz
import matplotlib.pyplot as plt

# read the data and set "datetime" as the index
bikes = pd.read_csv('../datasets/bikeshare.csv', index_col='datetime', parse_dates=True)

# "count" is a method, so it's best to rename that column
bikes.rename(columns={'count':'total'}, inplace=True)

# create "hour" as its own feature
bikes['hour'] = bikes.index.hour

bikes.head()

bikes.tail()

# - **hour** ranges from 0 (midnight) through 23 (11pm)
# - **workingday** is either 0 (weekend or holiday) or 1 (non-holiday weekday)

# # Exercise 7.1
#
# Run these two `groupby` statements and figure out what they tell you about the data.

# mean rentals for each value of "workingday"
bikes.groupby('workingday').total.mean()

# La cantidad de bicicletas promedio por día laboral y no laboral es similar, encontrando que es levemente superior la cantidad para el caso de días laborales.

bikes.groupby('hour').total.mean()

# mean rentals for each value of "hour"
bikes.groupby('hour').total.mean().plot(figsize=(15,5))

# Se identifican en general dos picos a largo del día, alrededor de las 8:00AM y las 5:00PM, lo cual podría obedecer a que en horas, la población utiliza la bicicleta para desplazasrse a trabajar o estudiar; es probable que a estas horas las personas también usen la bicicleta para ejercitarse. También, en la madrugada se observa la menor cantidad de bicicletas. A media mañana (alrededor de las 10:00 AM) encontramos un valle en los valores registrados.
#

# # Exercise 7.2
#
# Run this plotting code, and make sure you understand the output. Then, separate this plot into two separate plots conditioned on "workingday". (In other words, one plot should display the hourly trend for "workingday=0", and the other should display the hourly trend for "workingday=1".)

# mean rentals for each value of "hour"
bikes.groupby('hour').total.mean().plot()

# Plot for workingday == 0 and workingday == 1

# hourly rental trend for "workingday=0"
bikes[bikes.workingday==0].groupby('hour').total.mean().plot()

# hourly rental trend for "workingday=1"
bikes[bikes.workingday==1].groupby('hour').total.mean().plot()

# combine the two plots
bikes.groupby(['hour', 'workingday']).total.mean().unstack().plot()

# Write about your findings

# En workingday = 1 se observa el comportamiento de ingreso y salida de personas a su trabajo y en workingday=0 se observa solo una tendencia a salir por ser fin de semana/holiday/etc

# # Exercise 7.3
#
# Fit a linear regression model to the entire dataset, using "total" as the response and "hour" and "workingday" as the only features. Then, print the coefficients and interpret them. What are the limitations of linear regression in this instance?

# +
# Create linear regression object
regr = LinearRegression()

Y = bikes["total"]
X = bikes[["hour","workingday"]]

regr.fit(X, Y)
print('Coefficients: \n', regr.coef_)


# -

# Por cada unidad adicional de hora, el número total de bicicletas se incrementa en 10.49 unidades, dejando todo lo demás constante.
# Si el día no es holiday/weekend, el número total de bicicletas se incrementa en 4.07 unidades, dejando todo lo demás constante.
# Limitantes: no captura picos y valles (p.e. horas valle al final del día). No reconoce el comportamiento estacional/cíclico de la variable.

# # Exercice 7.4
#
# Create a Decision Tree to forecast "total" by manually iterating over the features "hour" and "workingday". The algorithm must at least have 6 end nodes.

def gini(y):
    if y.shape[0] == 0:
        return 0
    else:
        return 1 - (y.mean()**2 + (1 - y.mean())**2)


def gini_impurity(X_col, y, split):
    "Calculate the gain of an split k on feature j"
    
    filter_l = X_col < split
    y_l = y.loc[filter_l]
    y_r = y.loc[~filter_l]
    
    n_l = y_l.shape[0]
    n_r = y_r.shape[0]
    
    gini_y = gini(y)
    gini_l = gini(y_l)
    gini_r = gini(y_r)
    
    gini_impurity_ = gini_y - (n_l / (n_l + n_r) * gini_l + n_r / (n_l + n_r) * gini_r)
    
    return gini_impurity_


def best_split(X, y, num_pct=10):
    
    features = range(X.shape[1])
    
    best_split = [0, 0, 0]  # j, split, gain
    
    # For all features
    for j in features:
        
        splits = np.percentile(X.iloc[:, j], np.arange(0, 100, 100.0 / (num_pct+1)).tolist())
        splits = np.unique(splits)[1:]
        
        # For all splits
        for split in splits:
            gain = gini_impurity(X.iloc[:, j], y, split)
                        
            if gain > best_split[2]:
                best_split = [j, split, gain]
    
    return best_split


def tree_grow(X, y, level=0, min_gain=0.001, max_depth=None, num_pct=10):
    
    # If only one observation
    if X.shape[0] == 1:
        tree = dict(y_pred=y.iloc[:1].values[0], y_prob=0.5, level=level, split=-1, n_samples=1, gain=0)
        return tree
    
    # Calculate the best split
    j, split, gain = best_split(X, y, num_pct)
    
    # save tree and estimate prediction
    y_pred = int(y.mean() >= 0.5) 
    y_prob = (y.sum() + 1.0) / (y.shape[0] + 2.0)  # Laplace correction
    
    tree = dict(y_pred=y_pred, y_prob=y_prob, level=level, split=-1, n_samples=X.shape[0], gain=gain)
    
    # Check stooping criteria
    if gain < min_gain:
        return tree
    if max_depth is not None:
        if level >= max_depth:
            return tree   
    
    # No stooping criteria was meet, then continue to create the partition
    filter_l = X.iloc[:, j] < split
    X_l, y_l = X.loc[filter_l], y.loc[filter_l]
    X_r, y_r = X.loc[~filter_l], y.loc[~filter_l]
    tree['split'] = [j, split]

    # Next iteration to each split
    
    tree['sl'] = tree_grow(X_l, y_l, level + 1, min_gain=min_gain, max_depth=max_depth, num_pct=num_pct)
    tree['sr'] = tree_grow(X_r, y_r, level + 1, min_gain=min_gain, max_depth=max_depth, num_pct=num_pct)
    
    return tree


treebikes = tree_grow(X, Y, level=0, min_gain=0.001, max_depth=6, num_pct=10)

treebikes


def tree_predict(X, tree, proba=False):
    
    predicted = np.ones(X.shape[0])

    # Check if final node
    if tree['split'] == -1:
        if not proba:
            predicted = predicted * tree['y_pred']
        else:
            predicted = predicted * tree['y_prob']
            
    else:
        
        j, split = tree['split']
        filter_l = (X.iloc[:, j] < split)
        X_l = X.loc[filter_l]
        X_r = X.loc[~filter_l]

        if X_l.shape[0] == 0:  # If left node is empty only continue with right
            predicted[~filter_l] = tree_predict(X_r, tree['sr'], proba)
        elif X_r.shape[0] == 0:  # If right node is empty only continue with left
            predicted[filter_l] = tree_predict(X_l, tree['sl'], proba)
        else:
            predicted[filter_l] = tree_predict(X_l, tree['sl'], proba)
            predicted[~filter_l] = tree_predict(X_r, tree['sr'], proba)

    return predicted  


# Pte concluir Accuracy

# # Exercise 7.5
#
# Train a Decision Tree using scikit-learn. Comment about the performance of the models.

# +
import warnings
warnings.filterwarnings("ignore")

# list of values to try for max_depth
max_depth_range = range(1, 21)

# list to store the average RMSE for each value of max_depth
accuracy_scores = []

# use 10-fold cross-validation with each value of max_depth
from sklearn.model_selection import cross_val_score
from sklearn.tree import DecisionTreeClassifier

for depth in max_depth_range:
    clf = DecisionTreeClassifier(max_depth=depth, random_state=1)
    accuracy_scores.append(cross_val_score(clf, X, Y, cv=10, scoring='accuracy').mean())
# -

# plot max_depth (x-axis) versus RMSE (y-axis)
plt.plot(max_depth_range, accuracy_scores)
plt.xlabel('max_depth')
plt.ylabel('Accuracy')

# show the best accuracy and the corresponding max_depth
sorted(zip(accuracy_scores, max_depth_range))[::-1][0]

# max_depth=11 was best, so fit a tree using that parameter
clf = DecisionTreeClassifier(max_depth=11, random_state=1)
clf.fit(X, Y)

# # Part 2 - Bagging

# ## Mashable news stories analysis
#
# Predicting if a news story is going to be popular

df = pd.read_csv('../datasets/mashable.csv', index_col=0)
df.head()

df.shape

X = df.drop(['url', 'Popular'], axis=1)
y = df['Popular']

y.mean()

# train/test split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

y_train = pd.DataFrame(y_train)
y_test = pd.DataFrame(y_test)

# # Exercise 7.6
#
# Estimate a Decision Tree Classifier and a Logistic Regression
#
# Evaluate using the following metrics:
# * Accuracy
# * F1-Score

# +
from sklearn.linear_model import LogisticRegression
from sklearn import metrics

# Regression object
regr = LogisticRegression()
X = df.drop(['url', 'Popular'], axis=1)
y = df['Popular']

regr.fit(X, y)
y_pred = regr.predict(X_test)
print(metrics.f1_score(y_pred, y_test), metrics.accuracy_score(y_pred, y_test))

# +
from sklearn.tree import DecisionTreeClassifier

# Decision Tree
tree = DecisionTreeClassifier(max_depth=2)
X = df.drop(['url', 'Popular'], axis=1)
y = df['Popular']
tree.fit(X, y)
y_pred = tree.predict(X_test)
print(metrics.f1_score(y_pred, y_test), metrics.accuracy_score(y_pred, y_test))
# -

# # Exercise 7.7
#
# Estimate 300 bagged samples
#
# Estimate the following set of classifiers:
#
# * 100 Decision Trees where max_depth=None
# * 100 Decision Trees where max_depth=2
# * 100 Logistic Regressions

# +
n_estimators = 300
# set a seed for reproducibility
np.random.seed(123)

n_samples = X_train.shape[0]

samples = [np.random.choice(a=n_samples, size=n_samples, replace=True) for _ in range(n_estimators)]

# +
from sklearn.tree import DecisionTreeClassifier

np.random.seed(123) 
seeds = np.random.randint(1, 10000, size=n_estimators)
trees = {}
for i in range(n_estimators):
    if i < 100:
        trees[i] = DecisionTreeClassifier(max_features="sqrt", max_depth=None, random_state=seeds[i])
        trees[i].fit(X_train.iloc[samples[i]], y_train.iloc[samples[i]])
    if i >=200 :
        trees[i] = LogisticRegression()
        trees[i].fit(X_train.iloc[samples[i]], y_train.iloc[samples[i]])
    if i >=100 and i <= 199:
        trees[i] = DecisionTreeClassifier(max_features="sqrt", max_depth=2, random_state=seeds[i])
        trees[i].fit(X_train.iloc[samples[i]], y_train.iloc[samples[i]])
# -

y_pred_df = pd.DataFrame(index=X_test.index, columns=list(range(n_estimators)))
for i in range(n_estimators):
    y_pred_df.iloc[:, i] = trees[i].predict(X_test)
y_pred_df.head()

# # Exercise 7.8
#
# Ensemble using majority voting
#
# Evaluate using the following metrics:
# * Accuracy
# * F1-Score

y_pred_df.sum(axis=1)[:10]

y_pred_f= (y_pred_df.sum(axis=1) >= (n_estimators / 2)).astype(np.int)
from sklearn import metrics
print(metrics.f1_score(y_test, y_pred_f), metrics.accuracy_score(y_test, y_pred_f))

# # Exercise 7.9
#
# Estimate te probability as %models that predict positive
#
# Modify the probability threshold and select the one that maximizes the F1-Score

from sklearn import metrics
from sklearn.metrics import plot_precision_recall_curve
from sklearn.metrics import precision_recall_curve


y_pred_df.sum(axis=1)[:10]

y_pred_df = pd.DataFrame(index=X_test.index, columns=list(range(n_estimators)))
for i in range(n_estimators):
    y_pred_df.iloc[:, i] = trees[i].predict(X_test)

y_pred_f={}
for i in range(n_estimators):
    y_pred_df.iloc[:, i] = trees[i].predict(X_test)
    y_pred_f1 = trees[i].predict_proba(X_test)[:,1]
    y_pred_f[i] = y_pred_f1
y_pred_df.head()

# +
max_f1 = 0
max_tre = 0
for i in range(n_estimators):
    precision, recall, thresholds = precision_recall_curve(y_test, y_pred_f[i])
    f1_scores = 2 * recall *precision/(recall+precision)

    a = np.max(f1_scores)
    if max_f1 < np.max(f1_scores):
        max_f1 = np.max(f1_scores)
        max_tre =thresholds[np.argmax(f1_scores)]
        
print('Mejor Thresholds:',max_tre) 
print('Mejor F1:',max_f1)
# -

y_pred_2 = (y_pred_f1 >=0.1118499).astype(np.int)
metrics.f1_score(y_pred_2, y_test), metrics.accuracy_score(y_pred_2, y_test)

# # Exercise 7.10
#
# Ensemble using weighted voting using the oob_error
#
# Evaluate using the following metrics:
# * Accuracy
# * F1-Score

samples_oob = []
# show the "out-of-bag" observations for each sample
for sample in samples:
    samples_oob.append(sorted(set(range(n_samples)) - set(sample)))

# +
errors = np.zeros(n_estimators)

for i in range(n_estimators):
    y_pred_ = trees[i].predict(X_train.iloc[samples_oob[i]])
    errors[i] = 1 - metrics.accuracy_score(y_train.iloc[samples_oob[i]], y_pred_)
# -

y_train.iloc[samples_oob[1]]

# +
# %matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

plt.scatter(range(n_estimators), errors)
plt.xlim([0, n_estimators])
plt.title('OOB error of each tree')
# -

alpha = (1 - errors) / (1 - errors).sum()

# +
weighted_sum_1 = ((y_pred_df) * alpha).sum(axis=1)

y_pred = (weighted_sum_1 >= 0.5).astype(np.int)

metrics.f1_score(y_pred, y_test), metrics.accuracy_score(y_pred, y_test)
# -

# # Exercise 7.11
#
# Estimate te probability of the weighted voting
#
# Modify the probability threshold and select the one that maximizes the F1-Score

from sklearn.ensemble import BaggingClassifier
clf = BaggingClassifier(base_estimator=DecisionTreeClassifier(), n_estimators=300, bootstrap=True,
                        random_state=42, n_jobs=-1, oob_score=True)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
metrics.f1_score(y_pred, y_test), metrics.accuracy_score(y_pred, y_test)


# +
errors = np.zeros(clf.n_estimators)
y_pred_all_ = np.zeros((X_test.shape[0], clf.n_estimators))

for i in range(clf.n_estimators):
    oob_sample = ~clf.estimators_samples_[i]
    y_pred_ = clf.estimators_[i].predict(X_train.values[oob_sample])
    errors[i] = metrics.accuracy_score(y_pred_, y_train.values[oob_sample])
    y_pred_all_[:, i] = clf.estimators_[i].predict(X_test)
    
alpha = (1 - errors) / (1 - errors).sum()
y_pred = (np.sum(y_pred_all_ * alpha, axis=1) >= 0.5).astype(np.int)
# -

metrics.f1_score(y_pred, y_test), metrics.accuracy_score(y_pred, y_test)

# # Exercise 7.12
#
# Estimate a logistic regression using as input the estimated classifiers
#
# Modify the probability threshold such that maximizes the F1-Score

# +
from sklearn.linear_model import LogisticRegressionCV
from sklearn import metrics

X_train_2 = pd.DataFrame(index=X_train.index, columns=list(range(n_estimators)))
for i in range(n_estimators):
    X_train_2[i] = trees[i].predict(X_train)


# +
lr = LogisticRegressionCV(cv = 5 )
lr.fit(X_train_2, y_train)

print(metrics.f1_score(y_pred, y_test), metrics.accuracy_score(y_pred, y_test))

# -

y_pred_f1 = lr.predict_proba(y_pred_df)[:,1]

y_pred_f1

from sklearn.metrics import precision_recall_curve
from sklearn.metrics import plot_precision_recall_curve
precision, recall, thresholds = precision_recall_curve(y_test, y_pred_f1)

# +
f1_scores = 2*recall*precision/(recall+precision)

print('Mejor Thresholds:',thresholds[np.argmax(f1_scores)]) 
print('Mejor F1:',np.max(f1_scores))
# -

y_pred_2 = (y_pred_f1 >= 0.3081).astype(np.int)
metrics.f1_score(y_pred_2, y_test), metrics.accuracy_score(y_pred_2, y_test)



StarcoderdataJupyterScriptsDedupFiltered