Series temporales y ensembles de clasificadores con Python

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

data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month', header=0)
data.head()

data.plot(figsize=(15, 5))

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

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

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'])

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')

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')

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'])

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)

# Supongamos que los datos sin heterocedasticidad ni tendencia se llaman train_log_trend
plot_pacf(train_log_trend)
plt.show()

# Supongamos que los datos sin heterocedasticidad ni tendencia se llaman train_log_trend
plot_acf(train_log_trend)
plt.show()

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()

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')

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')

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')

# 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')

# 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"])

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)

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'])

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)

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()

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')

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)

# 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)

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)

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)

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)

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)

# 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)

# 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))

# 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))