Analisi dati con Python: data preprocessing e time series

A

Questa piccola guida ha come scopo di dimostrare in maniera pratica e comprensibile, come manipolare una serie di dati storici per far emergere rilevanti insights ed infine, fare una previsione futura utilizzando tecniche di machine learning, in particolare l’algoritmo ARIMA (autoregressive integrated moving average) questo modello è utilizzato molto spesso per fare previsioni su dati dove è presente una stagionalità.

L’algoritmo ARIMA consiste in una variazione di un modello ARMA con l’aggiunta di una componente d per la stazionarità, questo modello genera una serie di dati futuri utlizzando i suoi stessi dati passati.

Per essere impiegato richiede tre parametri: pd e q.

Con il primo parametro p ci si riferisce al il numero di lags usati nel modello, ovvero con che periodicità (seasonality) il modello ha a che fare. Il secondo parametro d si riferisce al grado di differenziazione, ovvero il numero di volte nelle quali i dati hanno avuto valori sottratti in passato mentre q rappresenta l’errore in forma di combinazione di errori passati.

Per chi fosse interessato ad approfondire, una spiegazione sicuramente più esaustiva e completa può essere trovata qui.

In questo articolo ho utilizzato Pyhton 3.6.5 in un Jupyter Notebookconfigurato tramite il gestore di pacchetti Miniconda.

N.B.

La prima parte di questo ‘tutorial’ è tranquillamente realizzabile con qualunque foglio di calcolo, quindi: non ho scoperto l’acqua calda.

Il dataset:

I dati utilizzati prevengono da un e-commerce specializzato in print on demand, quindi i datasets sono stati anonimizzati per rispettare la privacy dei clienti e, diverse colonne contenenti PIIs (personally identifiable information) sono state omesse.

Iniziamo importando le librerie necessarie e i datasets, dopodiché le colonne irrilevanti verrano rimosse dal dataset. Utilizzando pandas i dati vengono trasformati in un formato chiamato DataFrame , ovvero una sorta di array di arrays.

# Importazione delle librerie necessarie
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('ggplot')


# Importo i dati leggendo due file .csv tramite pandas
customers_df = pd.read_csv('customers_export.csv')
orders_df = pd.DataFrame(pd.read_csv('orders_export.csv'))
cols_to_drop_cust = [0, 1, 3, 4, 5, 6, 7, 9, 10, 12, 16, 17, 18]
customers_df.drop(customers_df.columns[cols_to_drop_cust], axis=1, inplace=True)

eseguendo queste poche righe di codice e chiamando il comando customers_df.head() avremo accesso alle prime cinque righe:

Una cosa che ha catturato subito la mia attenzione è stato vedere come più della metà degli indirizzi e-mail contenuti in questo dataset abbia un numero che spesso corrisponde all’anno di nascita. Quindi ho deciso di scrivere un piccolo algoritmo per risalire a questo numero, ed, eventualmente, creare un grafico con la distribuzione dei clienti per età.

# Calcolo l'età in base al numero estratto negli indirizzi mail
cust_age = customers_df.iloc[:,0].str.extract('(\d+)', expand=False).dropna().astype(int).tolist()
for n, i in enumerate(cust_age):
     if 70 <= i <= 100:
         cust_age[n] = i + 1900
for x in reversed(cust_age):
    if x > 2010:
        cust_age.remove(x)
for x in reversed(cust_age):
    if x < 70:
        cust_age.remove(x)
for x in reversed(cust_age):
    if x < 1900:
        cust_age.remove(x)

Dopo aver fatto un pò di pulizia e trasformato la lista in un oggetto Pandas.Series ho potuto create questo grafico, includendo solo le prime dieci leve:

Il prossimo step è quello di analizzare da che provincia i clienti ordinano, per questo basta selezionare la colonna relativa alla regione e utilizzando il metodo value_counts() :

# Clienti per regione
province = pd.DataFrame(customers_df.iloc[:,1])

# Elimino le righe contenenti dati 'null'
province = province.dropna(axis = 0, how= 'all')

province = province['Province Code'].value_counts()

ottengo questo ‘fantastico’ pie chart:

Sembrerebbe che dalle provincie di Milano e Roma proviene una grossa fetta di acquisti

Un altro punto interessante è tirare fuori quali sono le taglie più vendute, questo potrebbe aiutare a gestire in maniera migliore il magazzino, ed in un futuro, permetterebbe di utilizzare tecniche di machine learning per ottimizzare lo stesso.

#Convert list to df, drop NaN and plot 
tshirt_sizes_df = pd.DataFrame(tshirt_sizes)
tshirt_sizes_df = tshirt_sizes_df.dropna()

plt.plot(pd.value_counts(tshirt_sizes_df.values.flatten()))
plt.xlabel('Taglia')
plt.ylabel('Numero di t-shirt')
plt.title('Numero di t-shirt ordinate per taglia')

Analisi Time-series:

In questa parte analizzerò il dataset basandomi sulle date degli ordini, il dataset contiene informazioni dal 2016–05–01 al 2018–07–23, quindi più di due anni. Esistono diversi modelli time series, e per valutare quale conviene utilizzare è necessario fare uso di un toolkit che permette di scomporre la nostra serie nelle sue due componenti:

  • Sistematiche: componenti contenenti trend consistenti ed utilizzabili nel nostro modello
  • Non sistematiche: componenti della serie che non possono essere direttamente modellate

Qualunque serie di tipo time-series è composta da tre componenti sistematiche: livello, trend, stagionalità; e di una non sistematica, chiamata rumore.

Prima di poter utilizzare il modello però, è necessario fare un pò di ordine:

# selezione delle colonne
created_at = pd.read_csv(r'C:\Users\roberto.sannazzaro\Desktop\sales_2016-05-01_2018-07-23.csv')

drops = ['Sale ID', 'Order name', 'Transaction type', 'Sale type', 'Sales channel', 'Shipping country','Shipping region', 'Product variant SKU', 'POS location', 'Billing country',
        'Billing region', 'Billing city', 'Net quantity', 'Amount (before discount and taxes)', 'Line item discount', 'Order discount', 
         'Taxes', 'Product type', 'Product vendor', 'Amount (after discounts before tax)', 'Shipping city']

created_at.drop(drops, axis=1, inplace=True)
created_at = created_at.groupby(created_at['Date'])['Amount (after discounts and taxes)'].sum().reset_index()
created_at = created_at.set_index('Date')
created_at.index = pd.to_datetime(created_at.index)

y = created_at['Amount (after discounts and taxes)'].resample('M').mean()

# Creo il grafico
y.plot(figsize=(15, 6))
plt.xlabel('Data')
plt.ylabel('Media vendite in EUR')

 

Proiettando la serie si può notare come il trend è lievemente in crescita, quindi, come spiegato in questa lecture una scomposizione di tipo multiplicative fa al caso nostro, nel caso di trend con un crescita / perdita scarsa, oppure con una componente stagionale poco presente si può optare per un modello di tipo additive.

Per poter scomporre e creare il modello utilizzerò una biblioteca di python chiamata statsmodel .

import statsmodels.api as sm
from pylab import rcParams
rcParams['figure.figsize'] = 20, 12
decomposition = sm.tsa.seasonal_decompose(y, model='multiplicative', filt=None, two_sided=True, freq=7) # freq=1
fig = decomposition.plot()
plt.show()

Grazie a questa scomposizione ora, la componente stagionale è ben visibile e si può osservare come il trend in crescita. Però, è importante sottolineare che in questo caso, essendo a conoscenza delle dinamiche dietro questa attività commerciale ed essendo a conoscenza del fatto che nell’ultimo anno ci sono stati grandi cambiamenti in fatto di marketing e magazzino, opterò per considerare questo trend come neutro, in quanto non voglio che il modello venga influenzato da questi fattori esterni.

Il trend può essere positivo o negativo, lineare o non lineare, quindi è molto importante conoscere il proprio dataset per poter capire quando è passato abbastanza tempo per identificare un nuovo trend.

In questo momento, visti i risultati della scomposizione si può dire che ha senso utilizzare un modello ARIMA, per questo, ora è necessario scegliere i parametri P, D, Q (seasonality, trend, noise) per il modello.

Uno dei sistemi più rapidi ed efficaci per stimare questi parametri consiste nell’effettuare una ricerca chiamata “grid search” attraverso diverse combinazioni di valori P, D, Q combinate ad un indice di performance.

Per questo utilizzerò il “Test di verifica delle informazioni di Akaike” (AIC).

L’ AIC, data una serie di modelli misura la qualità di ciascuno in confronto agli altri modelli generati dalle combinazione di pd e q. Il valore dell’ AIC permetterà di capire come ciascun modello prende i dati in ingresso, tenendo conto la complessità del modello, quindi i modelli che avranno un risultato (fit) migliore, utilizzando meno features riceveranno un punteggio migliore (più basso).

la libreria per Pyhton pyramid-arima permette di utlizzare questo sistema di ricerca in maniera veloce, inoltre genera un oggetto con il quale si può direttamente chiamare il metodo fit() .

from pyramid.arima import auto_arima
stepwise_model = auto_arima(y, start_p=1, start_q=1,
                           max_p=3, max_q=3, m=12,
                           start_P=0, seasonal=True,
                           d=1, D=1, trace=True,
                           error_action='ignore',  
                           suppress_warnings=True, 
                           stepwise=True)
print(stepwise_model.aic())

Secondo le stime del test di Aike la combinazione ideale di variabili ha un punteggio di 112.10 con una combinazione di parametri 1, 1, 0, 12 .

Per poter verificare l’efficacia del modello, statsmodel fornisce una metodo chiamato plot_diagnostics() :

Non è perfetto, in quanto la linea contrassegnata con KDE non segue esattamente una distribuzione normale N(0, 1), questo si può anche osservare nel grafico Q-Q: la distribuzione dei residui non si discosta molto, per questo motivo le previsioni avranno un intervallo di confidenza maggiore (meno precise).

Una delle ragioni principali per la quale il modello non è perfetto consiste nella mancanza di dati, in quanto il dataset contiene informazioni dal momento esatto nel quale l’attivita commerciale ha iniziato, e quindi in una situazione instabile e senza alcun trend individuabile nei primi 6–8 mesi.

Validazione del modello:

Per validare il modello, e capire meglio l’accuratezza della previsioni, comparerò le vendite previste con quelle effettive, incominciando dal 2017–12–31:

Da quanto è possible osservare il modello prevede in maniera buona i dati reali ed è in grado di catturare la stagionalità presente nei dati.

Da quanto è possible osservare il modello prevede in maniera buona i dati reali ed è in grado di catturare la stagionalità presente nei dati.

Per avere una visione più precisa e misurabile è necessario calcolare il coefficiente R square, ovvero il quadrato della differenza tra valori stimati e quelli reali, più il valore di R square si avvicina a cento, migliori saranno le previsioni, ma è importante precisare che un valore basso di R square non sempre si associa ad un modello di bassa qualità, in quanto non esiste un range di valori standard nei quali il coefficiente R square si deve trovare, bensì dipende da caso a caso.

Allego due articoli che spigano in maniera più approfondita e completa il concetto di R square:

y_forecasted = pred.predicted_mean
y_truth = y['2017-01-01':]
mse = ((y_forecasted - y_truth) ** 2).mean()
print('R sqaure: {}'.format(round(mse, 2)))

Non male direi, anche se non è perfetto.

Come abbiamo potuto verificare il modello è in grado di catturare la stagionalità, ma più le previsioni si spingono nel futuro più l’intervallo di confidenza aumenterà. In questo esempio è presente solo una variabile: la media degli acquisti per settimana, mentre la variabile indipendente è descritta dalla data e dalla stagionalità ad essa associata.

In verde (speranza) il trend calcolato dal modello.

Per migliorare questo modello ed ottenere delle previsioni più precise, una buona idea sarebbe, ad esempio, includere una variabile esterna come il meteo: questo aiuterebbe a comprendere se c’è una correlazione tra meteo e acquisti e quindi, aiuterebbe (eventualmente) ad allocare il marketing budgetin maniera più efficace.

Complimenti per essere arrivato fin qui! Se sei interessato ad altri miei progetti puoi controllare la mia pagina oppure seguirmi su linkedin!

A proposito di me

roberto sannazzaro

I nostri Partner

Gli articoli più letti

Articoli recenti

Commenti recenti