In [ ]:
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

Carga y limpieza de datos

In [ ]:
df_start = pd.read_excel('../Data/Graficas DustTrack Honeywell 5sept comp.xlsx')
df_start.head()
In [ ]:
# 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))
In [ ]:
# Casting a variables numericas

df_start[['Valor DustTrack', 'Valor Honeywell']] = df_start[['Valor DustTrack', 'Valor Honeywell']].apply(pd.to_numeric)
In [ ]:
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')

Visualización inicial de los datos y exploraciones univariadas

In [ ]:
plot_time_series(df=df_start, variables=['Valor DustTrack', 'Valor Honeywell'], filename='Comparison_sensors')

Boxplot de ambas mediciones

In [ ]:
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
    
In [ ]:
build_boxplot(df=df_start, variables='Valor DustTrack', filename='Boxplot_DustTrack')
In [ ]:
build_boxplot(df=df_start, variables='Valor Honeywell', filename='Boxplot_Honeywell')
In [ ]:
build_boxplot(df=df_start, variables=['Valor DustTrack', 'Valor Honeywell'], filename='Boxplot_Comparison1')

Remoción de outliers

In [ ]:
df = df_start.copy()
df = df[np.logical_not(np.logical_or(df['Valor Honeywell'] > 90, df['Valor DustTrack'] > 492))]

Suavizamiento de serie con medias móviles

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

In [ ]:
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:]
In [ ]:
## Reescalamiento de los datos con la función Logaritmo
In [ ]:
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()
In [ ]:
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

In [ ]:
plot_time_series(df=movil_mean, variables=['Log10_DustTrack', 'Log10_Honeywell'], filename='log10_serie_movil')
In [ ]:
build_boxplot(df=movil_mean, variables=['Log10_DustTrack', 'Log10_Honeywell'], filename='Bloxplot_log10')

Partición de la base de datos en traint test

In [ ]:
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()

Escalamiento de los datos

In [ ]:
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()
In [ ]:
build_boxplot(df=df_train, variables=['scaled_DustTrack', 'scaled_Honeywell'], filename='Bloxplot_log10')

Correlación lineal de Pearson

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.

In [ ]:
# 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]
In [ ]:
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.

Construcción de modelos

In [ ]:
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))

Párametros de los modelos

In [ ]:
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]))

Append de predicciones a las series suavizadas

In [ ]:
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()

Append de predicciones a las series sin suaizar

In [ ]:
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()

Plots del valor real vs la predicción del modelo (sobre serie suavizada / train)

In [ ]:
plot_time_series(df=df_train, 
                 variables=['pred_min_max_scaler', 'pred_log_scaler', 'Mm_DustTrack'], 
                 filename='Predictions_01')

Plots del valor real vs la predicción del modelo (sobre serie suavizada / test)

In [ ]:
plot_time_series(df=df_test, 
                 variables=['pred_min_max_scaler', 'pred_log_scaler', 'Mm_DustTrack'], 
                 filename='Predictions_02')
In [ ]:
plot_time_series(df=df_start, 
                 variables=['pred_min_max_scaler', 'pred_log_scaler', 'Valor DustTrack'], 
                 filename='Predictions_03')

Errores de la predicción vs el valor real

In [ ]:
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

In [ ]:
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
In [ ]:
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')
In [ ]:
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')
In [ ]:
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')
In [ ]:
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']))
In [ ]:
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

Ecuación final de modelamiento (transformaión Log10)

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

Conclusiones finales

  • Es posible hacer una reconstrucción no perfecta pero si bastante buena de los datos del sensor DustTrack con base en los datos del sensor Honeywell.
In [ ]: