Series temporales y combinación de clasificadores

Esta práctica está dividida en dos partes:

  • En el primer ejercicio veremos como descomponer y componer series temporales para realizar predicciones a futuro.
  • En el segundo ejercicio estudiaremos diferentes métodos de combinación de clasificadores.
In [262]:
import pickle

import scipy.stats
import numpy as np
import pandas as pd
import matplotlib as mpl
from sklearn import svm
from sklearn import ensemble
from matplotlib import pyplot as plt
from sklearn.decomposition import PCA
from sklearn.metrics import accuracy_score
from sklearn.preprocessing import OneHotEncoder
from sklearn.tree import DecisionTreeClassifier
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
!pip install statsmodels
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_pacf, plot_acf

%matplotlib inline
Defaulting to user installation because normal site-packages is not writeable
Requirement already satisfied: statsmodels in ./.local/lib/python3.5/site-packages (0.11.1)
Requirement already satisfied: numpy>=1.14 in /usr/local/lib/python3.5/dist-packages (from statsmodels) (1.14.3)
Requirement already satisfied: scipy>=1.0 in ./.local/lib/python3.5/site-packages (from statsmodels) (1.4.1)
Requirement already satisfied: pandas>=0.21 in ./.local/lib/python3.5/site-packages (from statsmodels) (0.24.2)
Requirement already satisfied: patsy>=0.5 in ./.local/lib/python3.5/site-packages (from statsmodels) (0.5.1)
Requirement already satisfied: python-dateutil>=2.5.0 in /usr/local/lib/python3.5/dist-packages (from pandas>=0.21->statsmodels) (2.6.1)
Requirement already satisfied: pytz>=2011k in /usr/local/lib/python3.5/dist-packages (from pandas>=0.21->statsmodels) (2017.2)
Requirement already satisfied: six in /usr/local/lib/python3.5/dist-packages (from patsy>=0.5->statsmodels) (1.11.0)
WARNING: You are using pip version 20.0.2; however, version 20.1.1 is available.
You should consider upgrading via the '/usr/bin/python3 -m pip install --upgrade pip' command.

1. Series temporales

En este primer ejercicio trabajaremos las series temporales. Para ello usaremos el dataset AirPassangers que contiene información del número de vuelos que se realizaron a lo largo de muchos años.

Empezaremos leyendo los datos y observando gráficamente su distribución. Como se puede apreciar es un claro caso de serie temporal, con heterocedasticidad, tendencia, periodo y ruido. A lo largo de este ejercicio trataremos cada uno de estos puntos.

In [263]:
data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month', header=0)
data.head()
Out[263]:
Passengers
Month
1949-01-01112
1949-02-01118
1949-03-01132
1949-04-01129
1949-05-01121
In [264]:
data.plot(figsize=(15, 5))
Out[264]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff985d21a20>

Antes de empezar a tratar las diferentes componenetes de una serie temporal, eliminaremos del dataset original un par de años de datos. Así cuando hagamos una predicción a futuro podremos comprobar si se ajusta a los datos reales.

In [265]:
TEST_SIZE = 24
train, test = data.iloc[:-TEST_SIZE], data.iloc[-TEST_SIZE:]
x_train, x_test = np.array(range(train.shape[0])), np.array(range(train.shape[0], data.shape[0]))
train.shape, x_train.shape, test.shape, x_test.shape
Out[265]:
((120, 1), (120,), (24, 1), (24,))
In [266]:
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train)
ax.plot(x_test, test)
Out[266]:
[<matplotlib.lines.Line2D at 0x7ff9ccd10fd0>]

1.1 Componentes de la serie temporal

1.1.a. Heterocedasticidad

Se dice que una serie temporal tiene heterocedasticidad cuando la variancia varía con el tiempo (https://es.wikipedia.org/wiki/Heterocedasticidad). En nuestro caso, observamos que tenemos heterocedasticidad, ya que la amplitud de onda varía con el tiempo. En este primer apartado debéis eliminar la heterocedasticidad de la serie temporal. Es decir, que la diferencia entre el mínimo y el máximo de la estacionalidad (anual) sea más o menos la misma a lo largo del tiempo.

Implementación: Transformación de los datos para eliminar la heterocedasticidad. Para hacerlo debéis transformar los datos aplicando la función que consideréis que elimina mejor la heterocedasticidad en el dataset, provad con las tres funciones siguientes: - exponencial - logarítmica - raiz cuadrada y decidid cual de las tres funciona mejor para eliminar la heterocedasticidad. Estas funciones ya están implementadas en numpy. Mostrad gráficamente el dataset transformado.
In [267]:
train_log = np.log(train)
test_log = np.log(test)

print('logarítmica')
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_log)
ax.plot(x_test, test_log, label='test')
ax.set_title('logarítmica')
plt.legend([l1, l2],["train", 'test'])
logarítmica
Out[267]:
<matplotlib.legend.Legend at 0x7ff9ccb0edd8>
In [268]:
print('exponencial')
train_exp = np.exp(train)
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_exp)
ax.set_title('exponencial')
exponencial
Out[268]:
<matplotlib.text.Text at 0x7ff9cca66278>
In [269]:
print('raiz cuadrada')
train_sqrt = np.sqrt(train)
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_sqrt)
ax.set_title('raiz cuadrada')
raiz cuadrada
Out[269]:
<matplotlib.text.Text at 0x7ff9cca49e48>

Con el resultado obtenido, considero que la mejor opción para eliminarla heterocedasticidad en el dataset es la función logaritmica.

1.1.b. Tendencia

La tendencia es el comportamiento que tienen los datos a largo plazo (https://miro.medium.com/max/1872/1*rDQL2fAp_X_dgAHNZuwRfw.png). En nuestra serie temporal tenemos una tendencia lineal creciente. En este apartado debéis eliminar la tendencia, quedando una serie temporal con tendencia constante.

Implementación: Eliminar la tendencia de los datos. Como observando la serie podemos apreciar que tenemos una tendencia lineal, podemos ajustar una regresión lineal (usando scikit-learn) y sustraerla a los datos originales (sin heterocedasticidad). Graficar los datos sin tendencia.
In [270]:
model = LinearRegression()
model.fit(x_train.reshape(-1, 1), train_log.values)
train_log_trend = model.predict(x_train.reshape(-1, 1))



fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_log)
ax.plot(x_test, test_log)
ax.plot(x_train, train_log_trend)
plt.legend([l1, l2],["train", 'test'])

train_log_trend = (train_log - model.predict(x_train.reshape(-1, 1))).values.squeeze()


test_log_trend = (test_log - model.predict(x_test.reshape(-1, 1))).values.squeeze()


fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_log_trend)
ax.plot(x_test, test_log_trend)

plt.legend([l1, l2],["train", 'test'])
Out[270]:
<matplotlib.legend.Legend at 0x7ff99305ecc0>

1.1.c. Estacionalidad

Definimos estacionalidad como la variación cíclica que se produce en los datos (https://es.wikipedia.org/wiki/Estacionalidad). En este apartado se debe encontrar y eliminar la estacionalidad.

1.1.c.a. Encontrar el periodo de la estacionalidad

Implementación: En primer lugar tenéis que encontrar el ciclo, es decir, cada cuando se repiten los datos. Para encontrarlo podéis usar la autocorrelación (numpy te permite obtener los coeficientes de correlación). La idea es que se tiene que comparar: - La serie original con la serie empezando en el segundo punto (es decir, serie[1:]) - La serie original con la serie empezando en el tercer punto (es decir, serie[2:]) - ... En el momento en que vuelva a empezar el ciclo, la serie será muy parecida a la serie original, y, por lo tanto, la correlación será muy elevada. El punto donde la correlación sea máxima, será el ciclo de la serie. Mostrar gráficamente los coeficientes de los 22 primeros valores de autocorrelación y determinar su valor máximo, esta será nuestra estacionalidad. NOTA: Si usais numpy para encontrar los coeficientes de correlación, las dos series que comparéis deben de tener la misma longitud. Para representar la serie original podéis eliminar los valores del final para que tenga la misma longitud que la serie "movida". Por ejemplo, en el primer caso podéis comparar serie[1:] con serie[:-1].
In [271]:
corr_coefficients = []

for i in range(1, 24):
    corr_coefficients.append(np.corrcoef(train_log_trend[i:], train_log_trend[:-i])[0,1])

fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(range(len(corr_coefficients)), corr_coefficients)
Out[271]:
[<matplotlib.lines.Line2D at 0x7ff9ccc6f8d0>]

1.1.c.b. Aplicar un modelo SARIMA

Para tratar la estacionalidad aplicaremos un modelo SARIMA. Las siglas corresponden a "stationality ARIMA", es decir, un modelo ARIMA con estacionalidad (la que acabamos de encontrar).

El modelo ARIMA nos va a permitir tratar el ruido que queda al eliminar la heterocedasticidad, tendencia y estacionalidad. Para hacerlo tiene en cuenta las siguientes componentes:

  • AR: auto-regressive, se denomina p. Tiene en cuenta la correlación con sus lags, es decir, mira si las observaciones pasadas afectan para calcular el siguiente punto.
  • I: integrated, se denomina d. Es el orden de diferenciación, en nuestro caso no es necesario, puesto que la serie ya no tiene heterocedasticidad ni tendencia.
  • MA: moving-average, se denomina q. Tiene en cuenta la correlación con los lags de los errores, es decir, una vez se ha aplicado el modelo, mira los errores del modelo versus los datos reales.

Para más información sobre los ARIMA podéis consultar este link: https://otexts.com/fpp2/arima.html

Los modelos SARIMA dependen de varios parámetros (p,d,q)(P,D,Q)s, donde los primeros (minúsculas) corresponden a AR, I, MA de la serie normal, y los segundos (mayúsculas) a AR, I, MA con estacionalidad.

Para encontrar estos parámetros tenemos que mirar los gráficos PACF (partial autocorrelation function) y ACF (autocorrelation function).

El gráfico PACF nos determina el parámetro p, es decir, la AR.

In [272]:
# Supongamos que los datos sin heterocedasticidad ni tendencia se llaman train_log_trend
plot_pacf(train_log_trend)
plt.show()

Para leer este gráfico, simplemente nos tenemos que fijar en los valores que salen fuera del intervalo de confianza (zona azul).

Nada mas empezar hay dos valores fuera del intervalo. De todos modos el primero no se debe tener en cuenta, puesto que mira la correlación de un valor consigo mismo, y éste siempre será 1. Si no tenemos en cuenta este primer valor, solo hay un valor fuera del intervalo de confianza, con lo cual p = 1.

Cuando se repite el ciclo, es decir, a partir del valor 11, vuelven a haber dos valores fuera del intervalo, con lo cual P = 2.

Veamos ahora el gráfico ACF, este determinará el valor de q.

In [273]:
# Supongamos que los datos sin heterocedasticidad ni tendencia se llaman train_log_trend
plot_acf(train_log_trend)
plt.show()

Este gráfico se lee igual que el anterior. Nada más empezar hay un valor fuera del intervalo, con lo cual q = 1. Cuando se repite el ciclo, es decir, a partir del valor 11, hay tres valores fuera del intervalo, con lo cual Q = 3.

Como los datos no tienen ni tendencia ni heterocedasticidad, d = D = 0.

Como hemos visto en el apartado anterior, el ciclo es 12, con lo cual s = 12.

Implementación: Aplica un modelo SARIMA a los datos sin heterocedasticidad ni tendencia. Puedes usar SARIMAX (de statsmodels.tsa.statespace.sarimax) con los parámetros que acabamos de ver. Mostrar gráficamente el resultado del modelo junto con la serie original para comparar si se ajusta bien.
In [274]:
import warnings
warnings.filterwarnings("ignore") # specify to ignore warning messages
import statsmodels.api as sm
p=1
q=1
d=0
P=2
D=0
Q=3
m=12

#SARIMA(p,d,q)(P,D,Q)m

mod = sm.tsa.statespace.SARIMAX(train_log_trend,
                                order=(p,d,q),
                                seasonal_order=(P,D,Q, m))
results = mod.fit()
print('----Resultado del modelo----')
print(results.summary())

results.plot_diagnostics(figsize=(15, 12))

plt.show()
----Resultado del modelo----
                                         SARIMAX Results
==================================================================================================
Dep. Variable:                                          y   No. Observations:                  120
Model:             SARIMAX(1, 0, 1)x(2, 0, [1, 2, 3], 12)   Log Likelihood                 211.429
Date:                                    Sun, 07 Jun 2020   AIC                           -406.857
Time:                                            10:42:50   BIC                           -384.557
Sample:                                                 0   HQIC                          -397.801
                                                    - 120
Covariance Type:                                      opg
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.8580      0.078     11.035      0.000       0.706       1.010
ma.L1         -0.1886      0.122     -1.552      0.121      -0.427       0.050
ar.S.L12       0.1227      1.369      0.090      0.929      -2.560       2.805
ar.S.L24       0.8562      1.347      0.636      0.525      -1.783       3.496
ma.S.L12       0.3052      1.391      0.219      0.826      -2.420       3.031
ma.S.L24      -0.4774      0.787     -0.607      0.544      -2.019       1.065
ma.S.L36      -0.0355      0.190     -0.187      0.852      -0.407       0.336
sigma2         0.0013      0.000      7.820      0.000       0.001       0.002
===================================================================================
Ljung-Box (Q):                       38.61   Jarque-Bera (JB):                 3.61
Prob(Q):                              0.53   Prob(JB):                         0.16
Heteroskedasticity (H):               0.51   Skew:                             0.06
Prob(H) (two-sided):                  0.04   Kurtosis:                         3.84
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
In [275]:
print('En este punto, únicamente utilizo la parte test para ver si el modelo es fiel a la realidad')
start = len(train)
end = len(train) + len(test) - 1

predictions = results.predict(start, end,
                             typ = 'levels')

predictions= results.predict(start =start, end=end,exog=x_test)
forecast_1= results.forecast(steps=len(test) , exog=x_test)



train_size=len(train)


act= pd.DataFrame(test_log_trend)


predictions=pd.DataFrame(predictions)
predictions.reset_index(drop=True, inplace=True)
predictions.index=pd.DataFrame(x_test).index
predictions['Actual'] = act
predictions.rename(columns={0:'Pred'}, inplace=True)


predictions['Actual'].plot(figsize=(20,8), legend=True, color='blue')
predictions['Pred'].plot(legend=True, color='red', figsize=(20,8),title='Predición sobre el test')
En este punto, únicamente utilizo la parte test para ver si el modelo es fiel a la realidad
Out[275]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff991628a20>
In [276]:
print('En este punto, realizo una predicción sobre todo el dataset de entrenamiento')


start = 0
end = len(train) - 1

#https://medium.com/datadriveninvestor/time-series-prediction-using-sarimax-a6604f258c56
predictions = results.predict(start, end,
                             typ = 'levels')

predictions= results.predict(start =start, end=end)



train_size=0


act= pd.DataFrame(train_log_trend)


predictions=pd.DataFrame(predictions)
predictions.reset_index(drop=True, inplace=True)
predictions['Actual'] = act
predictions.rename(columns={0:'Pred'}, inplace=True)


predictions['Actual'].plot(figsize=(20,8), legend=True, color='blue')
predictions['Pred'].plot(legend=True, color='red', figsize=(20,8), title='dataset de entrenamiento')
En este punto, realizo una predicción sobre todo el dataset de entrenamiento
Out[276]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff991468780>
In [277]:
print('Utilizo la prediccion con forecast sobre el test')

forecast_apple= pd.DataFrame(forecast_1)
forecast_apple.reset_index(drop=True, inplace=True)
forecast_apple.index=pd.DataFrame(x_test).index
act= pd.DataFrame(test_log_trend)
forecast_apple['Actual'] =act
forecast_apple.rename(columns={0:'Forecast'}, inplace=True)


forecast_apple['Forecast'].plot(legend=True)
forecast_apple['Actual'].plot(legend=True, title='Predicción sobre el test con la función forecast')
Utilizo la prediccion con forecast sobre el test
Out[277]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff9913e3a20>

1.2. Predicción

Vamos a predecir los dos proximos años y comparar la predicción con los datos reales. Seguiremos los siguientes pasos:

  • Generar la predicción con el SARIMA
  • Añadir tendencia
  • Añadir heterocedasticidad

1.2.a. Preddicción SARIMA

Utilizaremos el modelo SARIMA que hemos ajustado antes para predecir los dos próximos años.

Implementación: Genera dos años de datos mediante el modelo SARIMA ajustado en el apartado anterior. Grafica toda la serie (sin heterocedasticidad ni tendencia), diferenciando con colores diferentes la serie real de los dos años de predicción.
In [278]:
# Forecast into future
print('Dos años son, los datos del test más  además los dos años que serían 24 meses')
print('No tenia claro si eran los dos años del dataset o los dos años más pero sin comparar. Pongo los dos por si acaso. ')
start = len(train)
end = len(train) + len(test) + 24
predictions = results.predict(start, end,
                             typ = 'levels')



act= pd.DataFrame(test_log_trend)


predictions=pd.DataFrame(predictions)
predictions.reset_index(drop=True, inplace=True)
predictions.index=pd.DataFrame(predictions).index
predictions['Actual'] = act
predictions.rename(columns={0:'Pred'}, inplace=True)


predictions['Actual'].plot(figsize=(20,8), legend=True, color='blue')
predictions['Pred'].plot(legend=True, color='red', figsize=(20,8), title='Predicción del test más dos años')
Dos años son, los datos del test más  además los dos años que serían 24 meses
No tenia claro si eran los dos años del dataset o los dos años más pero sin comparar. Pongo los dos por si acaso.
Out[278]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff9913475f8>
In [279]:
# Forecast for the next 2 years 
forecast = results.predict(start = len(train_log_trend),
                          end = (len(train_log_trend)-1) + 2 * 12,
                          typ = 'levels')

fig, ax = plt.subplots(1, 1, figsize=(15, 5))
l1,=ax.plot(range(0,len(train_log_trend)), train_log_trend)
ax.legend("train")

l2,=ax.plot(range(len(train_log_trend),(len(train_log_trend)) + 2 * 12), forecast)
ax.legend("predicción")


l3,=ax.plot(range(len(train_log_trend),(len(train_log_trend)) + 2 * 12), test_log_trend)
ax.legend("test")

plt.legend([l1, l2, l3],["train", 'predicción_test', "test_real"])
Out[279]:
<matplotlib.legend.Legend at 0x7ff991536828>

1.2.b. Tendencia

Anteriormente hemos visto que los datos de la serie temporal tienen una tendencia lineal y la hemos calculado mediante una regresión lineal. Vamos a añadir esta tendencia a nuestra predicción.

Implementación: Añadir a la serie anterior, la tendencia encontrada en el apartado 1.1.b. Mostrar gráficamente toda la serie (sin heterocedasticidad), diferenciando con colores diferentes la serie real de los dos años de predicción.
In [281]:
print('Este es un analisis personal sobre la tendencia.  No se solicita en el enunciado')
model = LinearRegression()
model.fit(x_train.reshape(-1, 1), train_log_trend)
train_log_trend_preduc = model.predict(x_train.reshape(-1, 1))

fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train_log_trend)
ax.plot(x_train, train_log_trend_preduc)

x=range(len(train_log_trend),(len(train_log_trend)) + 2 * 12)
ax.plot(x, forecast)




model = LinearRegression()
model.fit( np.array(list(x)).reshape(-1, 1), forecast)
train_log_trend_fores = model.predict(np.array(list(x)).reshape(-1, 1))
ax.plot(x, train_log_trend_fores)
Este es un analisis personal sobre la tendencia.  No se solicita en el enunciado
Out[281]:
[<matplotlib.lines.Line2D at 0x7ff9cc9f59e8>]
In [247]:
model = LinearRegression()
model.fit(x_train.reshape(-1, 1), train_log.values)

#sumo el valor
prediccion_log_sin_trend = (pd.DataFrame(forecast) + model.predict(x_test.reshape(-1, 1))).values.squeeze()

fig, ax = plt.subplots(1, 1, figsize=(15, 5))

test_trend_line= model.predict(x_test.reshape(-1, 1))

l1,=ax.plot(x_test, prediccion_log_sin_trend)
l2,=ax.plot(x_test, test_trend_line)


train_log_trend = model.predict(x_train.reshape(-1, 1))


l3,=ax.plot(x_test, test_log)
l4,=ax.plot(x_train, train_log)
l5,=ax.plot(x_train, train_log_trend)


plt.legend([l1, l2, l3, l4, l5],["predicción_con_tendencia", 'linea de la predción', "muestra_real", "entrenamiento",'linea_entrenamiento'])
Out[247]:
<matplotlib.legend.Legend at 0x7ff985d26ef0>

1.2.c. Heterocedasticidad

En el primer apartado de esta práctica hemos visto que la serie temporal tiene heterocedasticidad y la hemos eliminado transformando los datos. En este apartado haremos la transfromación inversa para añadir heterocedasticidad a nuestra predicción.

Implementació: Añadir a la serie anterior, la heterocedasticidad aplicando la función inversa a la encontrada en el apartado 1.1.a. Mostrar gráficamente toda la serie, diferenciando con colores diferentes la serie real de los dos años de predicción.
In [248]:
pred_inv_log = np.exp(prediccion_log_sin_trend)


fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_test, pred_inv_log)
Out[248]:
[<matplotlib.lines.Line2D at 0x7ff9ccbf0a58>]
Implementación: Añade al gráfico anterior la serie real para los dos años de predicción.
In [249]:
fig, ax = plt.subplots(1, 1, figsize=(15, 5))
ax.plot(x_train, train)
ax.plot(x_test, test, label='Real')
ax.plot(x_test, pred_inv_log, label='Prediction')
ax.legend()
Out[249]:
<matplotlib.legend.Legend at 0x7ff986b81898>

2. Combinación de clasificadores

Para este ejercicio usaremos el dataset cars.csv. Este dataset contiene datos sobre pruebas hechas en un modelo de coche, en concreto se estudia si el coche derrapará o no en una curva en función de la velocidad que lleva. El dataset contiene las siguientes columnas:

  • speed: velocidad que lleva el coche
  • angle: ángulo de la curva
  • succeed: 1 si el coche coge la curva sin problemas, 0 si el coche derrapa

A la primera parte de este ejercicio veremos la combinación de clasificadores en paralelo mediante las tecnicas de Bagging y Boosting.

La segunda parte pretende mejorar los resultados aplicando tecnicas de combinación secuencial de clasificadores: Stacking y Cascading.

Para empezar, vamos a visualizar el dataset.

In [250]:
cars = pd.read_csv('cars.csv')

fig, ax = plt.subplots(1, 1, figsize=(10, 5))
plt.scatter(cars['speed'], cars['angle'], c = cars['succeed'], cmap = 'bwr')
Out[250]:
<matplotlib.collections.PathCollection at 0x7ff986cdd438>

Para poder probar varios modelos, primero vamos a dividir el dataset entre train y test.

Para que todos obtengáis los mismos resultados y poder comentar dudas por el foro, fijaremos la seed para obtener los mismos datasets de train y test.

Como en este ejercicio trataremos stacking y cascading, y ambos se aplican sobre el conjunto de test, haremos un split del 50% para tener un poco más de base al aplicar estas dos técnicas.

In [251]:
X = cars.drop(['succeed'],1)
y = cars['succeed']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.5, random_state = 34)

Para visualizar gráficamente como funcionan los diferentes modelos, dibujaremos las fronteras de decisión. Lo que nos interesa es ver la forma que tienen (si es cuadriculada, suave, ...), y si hay zonas de indecisión (areas blancas). La idea es entender por qué el modelo que aplicamos tiene estas características.

Para dibujar las fronteras de decisión podéis usar la siguiente función.

In [252]:
# Creamis la meshgrid con los valores máximo y mínimo de 'x' e 'y'.
x_min, x_max = X['speed'].min()-0.1, X['speed'].max()+0.1
y_min, y_max = X['angle'].min()-0.1, X['angle'].max()+0.1

def plot_decision_boundaries(x, y, labels, model,
                             x_min=x_min,
                             x_max=x_max,
                             y_min=y_min,
                             y_max=y_max,
                             grid_step=1):

    xx, yy = np.meshgrid(np.arange(x_min, x_max, grid_step),
                         np.arange(y_min, y_max, grid_step))

    # Predecimos el clasificador con los valores de la meshgrid.
    Z = model.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:,1]

    # Hacemos un reshape para tener el formato correcto.
    Z = Z.reshape(xx.shape)

    # Seleccionamos una paleta de color.
    arr = plt.cm.coolwarm(np.arange(plt.cm.coolwarm.N))
    arr_hsv = mpl.colors.rgb_to_hsv(arr[:,0:3])
    arr_hsv[:,2] = arr_hsv[:,2] * 1.5
    arr_hsv[:,1] = arr_hsv[:,1] * .5
    arr_hsv = np.clip(arr_hsv, 0, 1)
    arr[:,0:3] = mpl.colors.hsv_to_rgb(arr_hsv)
    my_cmap = mpl.colors.ListedColormap(arr)

    # Hacemos el gráfico de las fronteras de decisión.
    fig, ax = plt.subplots(figsize=(7,7))
    plt.pcolormesh(xx, yy, Z, cmap=my_cmap)

    # Añadimos los puntos.
    ax.scatter(x, y, c=labels, cmap='coolwarm')
    ax.set_xlim(xx.min(), xx.max())
    ax.set_ylim(yy.min(), yy.max())
    ax.grid(False)

2.1. Combinación paralela de classificadores

2.1.1. Árboles de decisión

Para poder comparar el aumento de performance obtenido a medida que vamos aplicando técnicas nuevas, utilizaremos como baseline un simple árbol de decisión.

Implementación: Entrena un árbol de decisión sobre el conjunto de datos de train con profundidad máxima de 3 niveles (aplicaremos la misma restricción en las siguientes secciones). A continuación evalua sobre test y calcula su precisión aplicando validación cruzada con 5 conjuntos. Finalmente, dibuja la frontera de decisión. Sugerencia: usar el módulo *cross_val_score* de *sklearn*. Para aprender más sobre *cross validation* y sobre como usar estos módulos, os recomendamos los siguientes enlaces: - http://scikit-learn.org/stable/modules/cross_validation.html - http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_score.html
In [253]:
clf = DecisionTreeClassifier(max_depth=3)
cvscores = cross_val_score(clf, X_train, y_train, cv=5)
print("Precisión media obtenida con CV: {:.2f} +/- {:.2f} %".format(np.mean(cvscores)*100, np.std(cvscores)*100))
clf.fit(X_train, y_train)
plot_decision_boundaries(x=X_test['speed'], y=X_test['angle'], labels=y_test, model=clf)
Precisión media obtenida con CV: 91.34 +/- 1.38 %
Análisis: Analizad los resultados obtenidos y, en especial, la frontera de decisión.

Las fronteras de decisión están bastantemente definidas y observo que, apesar de algunos errores de clasificación se ha obtenido un modelo capaz de presentar la realidad. Las fronteras de decisión no són nada fluidas, eso es porque cada corte representa una división del árbol.

2.1.2.a. Bagging

La idea básica del bagging es utilizar el conjunto de entrenamiento original para generar centenares o miles de conjuntos similares usando muestreo con reemplazo. En este concepto está basado el algoritmo Random Forest, la combinación de varios árboles de decisión, cada uno entrenado con una realización diferente de los datos. La decisión final del clasificador combinado (la Random Forest) se toma por mayoría, dando el mismo peso a todas las decisiones parciales tomadas por los clasificadores base (los árboles).

Implementación: Entrena un random forest sobre el conjunto de datos de train con 20 árboles de decisión y profundidad máxima de 3 niveles. A continuación evalua sobre test y calcula su precisión aplicando validación cruzada con 5 conjuntos. Finalmente, dibuja la frontera de decisión. Sugerencia: usar el módulo *RandomForestClassifier* de *sklearn*. Per apender a usar este módulo os recomendamos el siguiente enlace: - http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html
In [254]:
from sklearn import ensemble

clf = ensemble.RandomForestClassifier(n_estimators=20, max_depth=3)

cvscores = cross_val_score(clf, X_train, y_train, cv=5)

print("Precisión media obtenida con CV: {:.2f} +/- {:.2f} %".format(np.mean(cvscores)*100, np.std(cvscores)*100))
clf.fit(X_train, y_train)


plot_decision_boundaries(x=X_test['speed'], y=X_test['angle'], labels=y_test, model=clf)
Precisión media obtenida con CV: 91.32 +/- 3.10 %
Análisis: Analizad los resultados obtenidos y, en especial, la frontera de decisión.

Mejoramos la precisión con respecto al módelo anterior. Al igual que el caso anterior, las fronteras de decisión no son nada fluidas, es decir, cada corte represent representa una división del árbol. Importante destacar la mejora (mínima) con respecto a un simple árbol de decisión.

2.1.3. Boosting

En el sistema de Boosting se combinan varios clasificadores débiles sequencialmente, y en cada uno de ellos se da más peso a los datos que han sido erróneamente clasificados en las combinaciones anteriores, para que se concentre así en los casos más difíciles de resolver.

Implementación: Entrena un gradient boosting sobre el conjunto de datos de train con 20 árboles de decisión y profundidad máxima de 3 niveles. A continuación evalua sobre test y calcula su precisión aplicando validación cruzada con 5 conjuntos. Finalmente, dibuja la frontera de decisión. Sugerencia: usar el módulo *GradientBoostingClassifier* de sklearn. Per apender a usar este módulo os recomendamos el siguiente enlace: - http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html
In [255]:
clf = ensemble.GradientBoostingClassifier(n_estimators=20, max_depth=3)

cvscores = cross_val_score(clf, X_train, y_train, cv=5)

print("Precisión media obtenida con CV: {:.2f} +/- {:.2f} %".format(np.mean(cvscores)*100, np.std(cvscores)*100))
clf.fit(X_train, y_train)


plot_decision_boundaries(x=X_test['speed'], y=X_test['angle'], labels=y_test, model=clf)
Precisión media obtenida con CV: 92.21 +/- 2.20 %
Análisis: Analizad los resultados obtenidos y, en especial, la frontera de decisión.

Mejoramos la precisión con respecto a los módelo anterior. Al igual que el caso anterior, las fronteras de decisión no son nada fluidas, es decir, cada corte represent representa una división del árbol. Importante destacar la mejora con respecto a un simple árbol de decisión. También destacar que tiene unas fronteras más marcadas que en un random forest.

2.2. Combinación secuencial de clasificadores base diferentes

Para poder hacer combinación secuencial de modelos, necesitamos tener varios modelos diferentes entrenados.

En nuestro caso, ya tenemos un árbol de decisión. Vamos a entrenar un par de modelos más.

2.2.1 KNN (k vecinos más próximos)

Implementación: Entrena un k-neighbors con 2 vecinos sobre el conjunto de datos de train. A continuación evalua sobre test y calcula su precisión aplicando validación cruzada con 5 conjuntos. Finalmente, dibuja la frontera de decisión.
In [256]:
clf = KNeighborsClassifier(n_neighbors=2)
clf.fit(X_train, y_train)

cvscores = cross_val_score(clf, X_train, y_train, cv=5)

print("Precisión media obtenida con CV: {:.2f} +/- {:.2f} %".format(np.mean(cvscores)*100, np.std(cvscores)*100))

plot_decision_boundaries(x=X_test['speed'], y=X_test['angle'], labels=y_test, model=clf)
Precisión media obtenida con CV: 89.18 +/- 4.13 %
Análisis: Analizad los resultados obtenidos y, en especial, la frontera de decisión.

Aunque la precisión del módelo es buena, no hemos obtenido una mejoria con respecto a los módelos anteriores. La frontera de decisión parece acurada, aunque hay algunas areas ("islas") que no clasifica bien.

2.2.2 SVM (Suport vector machines)

Implementació: Entrena un SVM con gamma = 0.07 sobre el conjunto de datos de train. A continuación evalua sobre test y calcula su precisión aplicando validación cruzada con 5 conjuntos. Finalmente, dibuja la frontera de decisión. NOTA: para que funcione la función de las fronteras de decisión se tiene que especificar "probability = True" al inicializar el modelo.
In [257]:
# Entrenem el classificador amb els paràmetres amb els que hem obtingut major precisió.
clf = svm.SVC( gamma=0.07, probability=True)
clf.fit(X_train, y_train)
cvscores = cross_val_score(clf, X_train, y_train, cv=5)

print("Precisión media obtenida con CV: {:.2f} +/- {:.2f} %".format(np.mean(cvscores)*100, np.std(cvscores)*100))

plot_decision_boundaries(x=X_test['speed'], y=X_test['angle'], labels=y_test, model=clf)
Precisión media obtenida con CV: 82.28 +/- 6.94 %
Análisis: Analizad los resultados obtenidos y, en especial, la frontera de decisión.

Aunque las fronteras de decisión son fluidas,se han observado demasiados islas y errores en la predicción. Este módelo es el que ha proporcionado un peor resultados.

2.2.1 Stacking

Un clasificador de stacking usa como atributos las predicciones hechas por otros clasificadores en lugar de los datos originales de entrada.

Implementación: Construye un clasificador de stacking usando un Gradient Boosting (con 20 árboles de decisión y profundidad máxima de 3 niveles) que use como atributos las predicciones hechas en el conjunto de test por los algoritmos: - árbol de decisión - knn - svm Calcula la precisión del modelo resultante con *cross-validation* en el conjunto de test (en este caso no tenemos conjunto de train, con lo cual se hace directamente cross-validation sobre test). Sugerencia: usar la función column_stack de numpy para juntar todas las predicciones. Per apender a usar esta función os recomendamos el siguiente enlace: - https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.column_stack.html
In [258]:
# Cargar las predicciones anteriores:
clf = DecisionTreeClassifier(max_depth=3)
clf.fit(X_train, y_train)
preds_tree = clf.predict(X_test)

clf = KNeighborsClassifier(n_neighbors=2, weights='uniform')
clf.fit(X_train, y_train)
preds_knn = clf.predict(X_test)

clf = svm.SVC( gamma=0.07, probability=True)
clf.fit(X_train, y_train)
preds_svm = clf.predict(X_test)

# Juntar las predicciones de los distintos clasificadores:
X_test_stacking = np.column_stack((preds_tree, preds_knn, preds_svm))
print("Dimensiones de la matriz con las predicciones de todos los clasificadores: {}".format(np.shape(X_test_stacking)))

# Transformar las variables categóricas (son todas) con OneHotEncoder:
enc = OneHotEncoder(n_values='auto')
enc.fit(X_test_stacking)
X_test_stacking = enc.transform(X_test_stacking).toarray()
print("Dimensiones de la matriz para entrenar el clasificador de stacking: {}".format(np.shape(X_test_stacking)))

# Calcular la precisión de un RandomForestClassifier con estas variables usando CV:
clf = ensemble.GradientBoostingClassifier(n_estimators=20, max_depth=3)
cvscores = cross_val_score(clf, X_test_stacking, y_test, cv=5)
print("Precisión media obtenida con CV: {:.2f} +/- {:.2f} %".format(np.mean(cvscores)*100, np.std(cvscores)*100))
Dimensiones de la matriz con las predicciones de todos los clasificadores: (232, 3)
Dimensiones de la matriz para entrenar el clasificador de stacking: (232, 6)
Precisión media obtenida con CV: 91.79 +/- 3.75 %
Análisis: ¿Has conseguido mejorar la precisión gracias al stacking? Comenta los resultados.

La precisión obtenida gracias al stacking es inferior a la obtenida con el mejor modelo anterior, gradient boosting (92.21). Otro detalle es la desviación de error, más elevada que la mejor opción.

2.2.2. Cascading

El caso de cascading es parecido al de stacking pero utilizando no solamente las predicciones parciales de los clasificadores base, sino también los datos originales.

Implementación: Construye un clasificador de cascading usando un Gradient Boosting (con 20 árboles de decisión y profundidad máxima de 3 niveles) que use como atributos las predicciones obtenidas con los modelos anteriores en el conjunto de test (igual que con el stacking), y también las variables originales. Calcula la precisión del modelo resultante con *cross-validation* en e conjunto de test. Sugerencia: Usa el mismo conjunto de datos que en el ejercicio anterior pero añade `X_test`.
In [259]:
# Juntar las predicciones de los distintos clasificadores y las variables originales:
X_test_cascading = np.column_stack((X_test, X_test_stacking))
print("Dimensiones de la matriz: {}".format(np.shape(X_test_cascading)))

# Calcular la precisión de un RandomForestClassifier con estas variables usando CV:
clf = ensemble.GradientBoostingClassifier(n_estimators=20, max_depth=3)
cvscores = cross_val_score(clf, X_test_cascading, y_test, cv=5)
print("Precisión media obtenida con CV: {:.2f} +/- {:.2f} %".format(np.mean(cvscores)*100, np.std(cvscores)*100))
Dimensiones de la matriz: (232, 8)
Precisión media obtenida con CV: 93.52 +/- 4.08 %
Análisis: ¿Has conseguido mejorar la precisión gracias al cascading? Comenta los resultados.

La precisión obtenida gracias al cascading es ligeramente superior a la obtenida con stacking y superior a la mejor opción hasta el momento. Las desviaciones de las precisiones medias con CV también ha sido superior.