import dateutil
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import plotly.plotly as py
import plotly.graph_objs as go
import plotly.figure_factory as ff
import statsmodels.api as sm
from scipy import stats
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn import linear_model
from sklearn.preprocessing import MinMaxScaler
df_start = pd.read_excel('../Data/Graficas DustTrack Honeywell 5sept comp.xlsx')
df_start.head()
# Eliminar nulos y headers repetidos
df_start = df_start[~df_start['Valor DustTrack'].isna()]
df_start = df_start[df_start.columns.values[:-1]]
sum(df_start['Valor DustTrack'] == 'Valor DustTrack')
df_start = df_start[np.logical_not(df_start['Valor DustTrack'] == 'Valor DustTrack')]
df_start.index = np.arange(len(df_start))
# Casting a variables numericas
df_start[['Valor DustTrack', 'Valor Honeywell']] = df_start[['Valor DustTrack', 'Valor Honeywell']].apply(pd.to_numeric)
def plot_time_series(df, variables, filename):
data = list()
for variable in variables:
data.append(go.Scatter(x = df.index,
y = df[variable],
name=variable))
layout = go.Layout(title='DustTrack vs Honeywell',
yaxis=dict(title='Pm 2.5 value'),
xaxis=dict(title='N° Observation'))
fig = go.Figure(data=data, layout=layout)
plot_url = py.iplot(fig, filename=filename)
return plot_url
# Login plotly
py.sign_in('FranciscoJavierOspinaSalazar', 'WLz5dfnwRjfo2cSwlRiL')
plot_time_series(df=df_start, variables=['Valor DustTrack', 'Valor Honeywell'], filename='Comparison_sensors')
def build_boxplot(df, variables, filename):
if isinstance(variables, str):
variables = [variables]
data = list()
for variable in variables:
data.append(go.Box(y = df[variable],
name=variable))
layout = go.Layout(title='Boxplot',
yaxis=dict(title='Pm 2.5 value'))
fig = go.Figure(data=data, layout=layout)
plot_url = py.iplot(fig, filename=filename)
return plot_url
build_boxplot(df=df_start, variables='Valor DustTrack', filename='Boxplot_DustTrack')
build_boxplot(df=df_start, variables='Valor Honeywell', filename='Boxplot_Honeywell')
build_boxplot(df=df_start, variables=['Valor DustTrack', 'Valor Honeywell'], filename='Boxplot_Comparison1')
df = df_start.copy()
df = df[np.logical_not(np.logical_or(df['Valor Honeywell'] > 90, df['Valor DustTrack'] > 492))]
La descripción de las médias móviles y su uso para suavizar picos en series de tiempo se encuentra en http://www.expansion.com/diccionario-economico/media-movil.html
n_movil_mean = 20
movil_mean = df[['Valor DustTrack', 'Valor Honeywell']].rolling(n_movil_mean).mean()
movil_mean.columns = ['Mm_DustTrack', 'Mm_Honeywell']
movil_mean = movil_mean.loc[n_movil_mean - 1:]
## Reescalamiento de los datos con la función Logaritmo
movil_mean = movil_mean.assign(Log10_DustTrack=np.log10(movil_mean['Mm_DustTrack']))
movil_mean = movil_mean.assign(Log10_Honeywell=np.log10(movil_mean['Mm_Honeywell']))
movil_mean.index = np.arange(len(movil_mean))
movil_mean.head()
plot_time_series(df=movil_mean, variables=['Mm_DustTrack', 'Mm_Honeywell'], filename='serie_movil')
Con el suavizamiento se suprimieron los picos que superaban 2000 en el eje y
plot_time_series(df=movil_mean, variables=['Log10_DustTrack', 'Log10_Honeywell'], filename='log10_serie_movil')
build_boxplot(df=movil_mean, variables=['Log10_DustTrack', 'Log10_Honeywell'], filename='Bloxplot_log10')
size_train = int(len(movil_mean) * .7)
index_train = np.random.choice(np.arange(len(movil_mean)), size=size_train, replace=False)
index_test = np.setdiff1d(np.arange(len(movil_mean)), index_train)
df_train = movil_mean.iloc[index_train]
df_test = movil_mean.iloc[index_test]
df_train = df_train.sort_index()
df_test = df_test.sort_index()
scaler_dt = MinMaxScaler()
scaler_dt.fit(df_train['Mm_DustTrack'].values.reshape(-1, 1))
scaled_DustTrack_train = scaler_dt.transform(df_train['Mm_DustTrack'].values.reshape(-1,1))
scaled_DustTrack_test = scaler_dt.transform(df_test['Mm_DustTrack'].values.reshape(-1,1))
scaler_hw = MinMaxScaler()
scaler_hw.fit(df_train['Mm_Honeywell'].values.reshape(-1, 1))
scaled_Honeywell_train = scaler_hw.transform(df_train['Mm_Honeywell'].values.reshape(-1,1))
scaled_Honeywell_test = scaler_hw.transform(df_test['Mm_Honeywell'].values.reshape(-1,1))
df_train = df_train.assign(scaled_DustTrack=scaled_DustTrack_train)
df_train = df_train.assign(scaled_Honeywell=scaled_Honeywell_train)
df_test = df_test.assign(scaled_DustTrack=scaled_DustTrack_test)
df_test = df_test.assign(scaled_Honeywell=scaled_Honeywell_test)
movil_mean = movil_mean.assign()
df_train.head()
build_boxplot(df=df_train, variables=['scaled_DustTrack', 'scaled_Honeywell'], filename='Bloxplot_log10')
https://es.wikipedia.org/wiki/Coeficiente_de_correlaci%C3%B3n_de_Pearson
Debido a que solo son dos variables si la correlación de pearson es alta se pueden hacer modelos que intenten generar la variable Valor DustTrack a partir de Valor Honeywell.
# Coeficinte de correlación sobre las variables transformadas por logaritmos de las
# series suavizadas por las media móviles
stats.pearsonr(df_train['Log10_DustTrack'], df_train['Log10_Honeywell'])[0]
stats.pearsonr(df_train['scaled_DustTrack'], df_train['scaled_Honeywell'])[0]
Tanto en la gráfica de la serie de tiempo como en la correlación, estas son las mejores variables hasta el momento para el modelo.
model_min_max_scaler = linear_model.LinearRegression(fit_intercept=True)
model_log_scaler = linear_model.LinearRegression(fit_intercept=True)
model_min_max_scaler.fit(X=df_train['scaled_Honeywell'].values.reshape(-1,1), y=df_train['scaled_DustTrack'].values.reshape(-1,1), )
model_log_scaler.fit(X=df_train['Log10_Honeywell'].values.reshape(-1,1), y=df_train['Log10_DustTrack'].values.reshape(-1,1))
print('Parameters model with log scaler. Coef = {}, Intercept = {}'.format(model_log_scaler.coef_[0][0], model_log_scaler.intercept_[0]))
print('Parameters model with min_max scaler. Coef = {}, Intercept = {}'.format(model_min_max_scaler.coef_[0][0], model_min_max_scaler.intercept_[0]))
pred_min_max_scaler_train = scaler_dt.inverse_transform(model_min_max_scaler.predict(df_train['scaled_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
pred_log_scaler_train = np.power(10, model_log_scaler.predict(df_train['Log10_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
df_train = df_train.assign(pred_min_max_scaler=pred_min_max_scaler_train)
df_train = df_train.assign(pred_log_scaler=pred_log_scaler_train)
pred_min_max_scaler_test= scaler_dt.inverse_transform(model_min_max_scaler.predict(df_test['scaled_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
pred_log_scaler_test = np.power(10, model_log_scaler.predict(df_test['Log10_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
df_test = df_test.assign(pred_min_max_scaler=pred_min_max_scaler_test)
df_test = df_test.assign(pred_log_scaler=pred_log_scaler_test)
df_test.head()
df_start = df_start.assign(pred_log_scaler=np.power(10, model_log_scaler.predict(np.log10(df_start['Valor Honeywell']).values.reshape(-1,1)).reshape(-1,1)))
scaled_Honeywell = scaler_hw.transform(df_start['Valor Honeywell'].values.reshape(-1,1))
df_start = df_start.assign(scaled_Honeywell=scaled_Honeywell)
pred_min_max_scaler = scaler_dt.inverse_transform(model_min_max_scaler.predict(df_start['scaled_Honeywell'].values.reshape(-1,1)).reshape(-1,1))
df_start = df_start.assign(pred_min_max_scaler=pred_min_max_scaler)
df_start.head()
plot_time_series(df=df_train,
variables=['pred_min_max_scaler', 'pred_log_scaler', 'Mm_DustTrack'],
filename='Predictions_01')
plot_time_series(df=df_test,
variables=['pred_min_max_scaler', 'pred_log_scaler', 'Mm_DustTrack'],
filename='Predictions_02')
plot_time_series(df=df_start,
variables=['pred_min_max_scaler', 'pred_log_scaler', 'Valor DustTrack'],
filename='Predictions_03')
dc_errors = dict()
for key, value in {'Train Data': df_train, 'Test Data': df_test, 'Original Data': df_start}.items():
true_val = 'Mm_DustTrack'
if key == 'Original Data':
true_val = 'Valor DustTrack'
dc_errors[key] = [mean_squared_error(value[true_val], value['pred_log_scaler']),
mean_squared_error(value[true_val], value['pred_min_max_scaler']),
mean_absolute_error(value[true_val], value['pred_log_scaler']),
mean_absolute_error(value[true_val], value['pred_min_max_scaler'])]
df_errors = pd.DataFrame(dc_errors)
df_errors.index = ['Log10 transform MSE', 'Min_max transform MSE', 'Log10 transform MAE', 'Min_max transform MAE']
df_errors[['Train Data', 'Test Data', 'Original Data']]
Nota: Los valores de MAE y MSE corresponden al error absoluto medio y error cuadrático medio
def plot_density(df, y_tue, y_pred, group_labels, title, filename):
hist_data = list()
for y in y_pred:
hist_data.append(df[y_tue] - df[y])
group_labels = group_labels
# Create distplot with custom bin_size
fig = ff.create_distplot(hist_data, group_labels, show_hist=False)
fig['layout'].update(title=title)
plot_url = py.iplot(fig, filename=filename)
return plot_url
plot_density(df=df_train,
y_tue='Mm_DustTrack',
y_pred=['pred_log_scaler', 'pred_min_max_scaler'],
group_labels=['Log scaler', 'Min_max scaler'],
title='Densidad de los residuales train',
filename='Residuals_density_01')
plot_density(df=df_test,
y_tue='Mm_DustTrack',
y_pred=['pred_log_scaler', 'pred_min_max_scaler'],
group_labels=['Log scaler', 'Min_max scaler'],
title='Densidad de los residuales test',
filename='Residuals_density_02')
plot_density(df=df_start,
y_tue='Valor DustTrack',
y_pred=['pred_log_scaler', 'pred_min_max_scaler'],
group_labels=['Log scaler', 'Min_max scaler'],
title='Densidad de los residuales test',
filename='Residuals_density_02')
df_train = df_train.assign(res_log_scaler=np.abs(df_train['Mm_DustTrack'] - df_train['pred_log_scaler']))
df_train = df_train.assign(res_min_max_scaler=np.abs(df_train['Mm_DustTrack'] - df_train['pred_min_max_scaler']))
df_test = df_test.assign(res_log_scaler=np.abs(df_test['Mm_DustTrack'] - df_test['pred_log_scaler']))
df_test = df_test.assign(res_min_max_scaler=np.abs(df_test['Mm_DustTrack'] - df_test['pred_min_max_scaler']))
df_start = df_start.assign(res_log_scaler=np.abs(df_start['Valor DustTrack'] - df_start['pred_log_scaler']))
df_start = df_start.assign(res_min_max_scaler=np.abs(df_start['Valor DustTrack'] - df_start['pred_min_max_scaler']))
step = 0.05
quantile_train = df_train[['res_log_scaler', 'res_min_max_scaler']].quantile(np.arange(0, 1 + step, step))
quantile_train.columns = [['res_log_scaler_train', 'res_min_max_scaler_train']]
quantile_test = df_test[['res_log_scaler', 'res_min_max_scaler']].quantile(np.arange(0, 1 + step, step))
quantile_test.columns = [['res_log_scaler_test', 'res_min_max_scaler_test']]
quantile_start = df_start[['res_log_scaler', 'res_min_max_scaler']].quantile(np.arange(0, 1 + step, step))
quantile_start.columns = [['res_log_scaler_original_data', 'res_min_max_scaler_original_data']]
df_quantile = pd.concat([quantile_train, quantile_test, quantile_start], axis=1)
df_quantile
La ecuación final que mejor se comporta para el modelamiento de los valores de Valor DustTrack a partir de los datos de Valor Honeywell es:
$$y = 10^{1.1675 * log_{10}(x) + 0.4173}$$Donde X es el valor de la variable Valor Honeywell y y el valor de la variable Valor DustTrack m