Breve introduzione al Dynamic Time Warping

Overview

1. Dynamic Time Warping
   1.1 Definizione
   1.2 Algoritmo
   1.3 Window Constraint
   1.4 Miglioramenti al DTW
2. DTW e Machine Learning

Image Credits

Il codice di questo articolo è disponibile su un notebook Google Colab.

Dynamic Time Warping

Definizione

Dynamic Time Warping (DTW) è un algoritmo utilizzato per misurare la similarità tra due o più serie temporali, andando a leggere su Wikipedia troviamo la seguente definizione:

“Il Dynamic time warping, o DTW, è un algoritmo che permette l’allineamento tra due sequenze, e che può portare ad una misura di distanza tra le due sequenze allineate. Tale algoritmo è particolarmente utile per trattare sequenze in cui singole componenti hanno caratteristiche che variano nel tempo, e per le quali la semplice espansione o compressione lineare delle due sequenze non porta risultati soddisfacenti.”

Che casi d’uso ci suggerisce la definizione? Prendiamo ad esempio il riconoscimento vocale, una persona potrebbe parlare più velocemente o più lentamente della voce campione preregistrata, il DTW non risente di questo problema a differenza di una distanza che venga applicata “punto per punto”.

Il DTW è molto utilizzato in ambito IoT dove gli stream di dati proveninenti da video, audio, sensori, possono essere trasformati in una sequenza lineare per essere analizzati con questa tecnica, in particolare è molto efficace con serie ad alta frequenza di campionamento, pensiamo ad esempio al mercato delle utilities o dell’automotive.

L’idea di base è calcolare la distanza tra due vettori/serie storiche indipendentemente dalla loro lunghezza e per farlo bisogna che i vettori rispettino le seguenti regole:

• ogni indice della prima sequenza deve poter essere confrontato con tutti gli indici delle altre e viceversa
• il primo indice della prima sequenza deve essere associato al primo indice dell’altra sequenza (ma non deve essere la sua unica corrispondenza)
• l’ultimo indice della prima sequenza deve corrispondere all’ultimo indice dell’altra sequenza (ma non deve essere la sua unica corrispondenza)
• La mappatura degli indici dalla prima sequenza agli indici dell’altra sequenza deve essere monotonicamente crescente. A grandi linee, indipendentemente dalla lunghezza delle serie l’indice iniziale e finale devono coincidere

Algoritmo

Per illustrare il DTW userò in parte un breve esempio con due sequenze di pochi elementi preso da riptutorial.

L’algoritmo di base è molto semplice, vedremo la differenza della similarità tra serie con DTW e distanza euclidea, partiamo creando gli array.

import numpy as np
serie_1 = np.array([1, 2, 3, 5, 5, 5, 6], ndmin = 2)
serie_2 = np.array([1, 1, 2, 2, 3, 5], ndmin =2)

Le due serie hanno lunghezza diversa quindi un confronto punto per punto con distanza euclidea (o altra distanza) non sarebbe del tutto corretto. Supponiamo che le due serie partano nello stesso periodo e che la serie_2 finisca prima, quindi aggiungiamo un elemento per calcolare la distanza (non rispetta le regole sugli indici del DTW).

\[ d\left( p,q\right) = \sqrt {\sum _{i=1}^{n} \left( q_{i}-p_{i}\right)^2 } \]

serie_2New = serie_2.copy()
serie_2New = np.append(serie_2New, 0)

distance = np.sqrt(np.sum((serie_1-serie_2New)**2))
print(distance)

La distanza calcolata è 7.14.

Passiamo ora al calcolo della distanza mediante DTW, per farlo ci servirà una matrice di dimensione n x m dove m ed n sono la lunghezza della serie più 1.

n = serie_1.shape[1]+ 1
m = serie_2.shape[1]+ 1
DTW = np.zeros([n,m])

Il primo elemento della matrice DTW[0,0] deve essere pari a 0, mentre prima riga e prima colonna pari ad infinito. Quanto fatto serve per il cold start dell’algoritmo.
Il cold start è il confronto con i periodi precedenti quando si è al primo passo dell’algoritmo, impostando una riga fittizia ad infinito non si incorre in situazioni anomale e si evita di implementare condizioni sullo start che appesantirebbero ulteriormente il calcolo.

ecco spiegato perchè la dimensione della serie più 1.

DTW[0,1:] = np.inf
DTW[1:,0] = np.inf
print(DTW)

Per ogni posizione della matrice verrà calcolata la distanza euclidea di ogni punto di una serie da tutti i punti delle altre partendo dalla posizione DTW[1,1]. Ad ogni distanza viene aggiunto il minimo tra:

• il valore precedente sulle righe: cancellazione
• il valore precedente sulle colonne: inserimento
• il valore precedente in diagonale: corrispondenza

for i in range(1,n):
    for j in range(1, m):
        dist = abs(serie_1[:,i-1] - serie_2[:,j-1]) 
        DTW[i,j] = dist + min([DTW[i-1,j], DTW[i,j-1], DTW[i-1,j-1]])

print(DTW)

L’algoritmo restituisce una matrice con tutte le posizioni valorizzate.
Arrivati a questo punto partendo da DTW[0,0] bisogna costruire un percorso che tocchi il valore minimo più vicino a quello nella posizione in cui ci troviamo muovendoci da sinistra verso destra in orizzontale/verticale/obliquo, questo è il Warping Path.

• uno spostamento orizzontale significa che serie_2 è accelerata durante questo intervallo.
• uno spostamento verticale significa che serie_2 è decelerata durante questo intervallo.
• una mossa diagonale significa che durante questo periodo le serie camminano di pari passo

Il valore finale del path è la distanza tra le serie.
Dalla matrice si nota che il calcolo non è stato influenzato dalla diversa dimensione delle sequenze e che il valore finale è molto più piccolo di quello calcolato punto per punto con la distanza euclidea.

Window Constraint

La complessità di calcolo del DTW è O(m * n) dove m e n rappresentano la lunghezza di ciascuna sequenza quindi se avessimo sequenze molto lunghe il confronto tra le due risulterebbe complesso, ed anche inutile, perchè alcuni confronti verrebbero effettuati su posizioni lontanissime.
Prendiamo l’esempio della serie storica, vorrebbe dire confrontare periodi molto distanti temporalmente che non hanno nessuna relazione ed influenza tra loro.
Per questo motivo viene inserita un parametro w, una finestra temporale su cui calcolare le differenze.

La finestra temporale è il valore minimo tra un intero definito dall’utente (io ho utilizzato 1) e la differenza in valore assoluto tra m ed n.

w = max([1, abs(n-m)])

Di seguito come cambia l’algoritmo per il calcolo della distanza.

for i in range(1,n):
    for j in range(max([1,i-w]), min([m, i+w])):
        dist = abs(serie_1[:,i-1] - serie_2[:,j-1]) 
        DTW[i,j] = dist + min([DTW[i-1,j], DTW[i,j-1], DTW[i-1,j-1]])

print(DTW)

Nell’esempio proposto non si noteranno differenze viste le sequenze molto corte e con una differenza di una sola posizione, ma per sequenze molto lunghe e con differenze di lunghezza evidenti l’impatto è tangibile.

Miglioramenti al DTW

Vista la complessità quadratica del DTW sono stati proposti alcuni migliorameti dell’algoritmo, il più conosciuto è probabilmente il FastDTW.
Il FastDTW cerca di migliorare le performance approssimando alcuni risultati, non sempre questa riduzione di tempi si traduce in una qualità del risultato accettabile.
Ci sono molti paper e package su questo algoritmo, vi lascio il link della libreria python fastdtw e il paper “FastDTW is approximate and Generally Slower than the Algorithm it Approximates” che spiega in quali casi questo algoritmo sia davvero più performante (mantenendo accettabile l’efficacia) del DTW classico.

DTW e Machine Learning

Tornando a quanto detto nell’introduzione, come può essere/viene usato il DTW nei problemi di Machine Learning?

• apprendimento supervisionato, come distanza per algoritmi KNN su Time Series
• apprendimento non supervisionato, come misura di similarità nel clustering di Time Series
• apprendimento semi-supervisionato, etichettare le Time Series utilizzando un algoritmo di clustering basato su DTW e successivamente classificare secondo queste assegnazioni

Di seguito un esempio di clustering basato sul DTW.

Ingrandendo il grafico finale si può notare come alcune serie dei cluster B e C sembrino molto più vicine al centroide A in alcuni periodi ma il DTW non ne risente e permette un raggruppamento migliore.

Lascio alcuni paper del 2020 sull’argomento per chi volesse approfondire:

• Time Series Clustering Model Based on DTW for Classifying Car Parks
• Intertemporal Similarity of Economic Time Series: An Application of Dynamic Time Warping
• Time series classification using local distance-based features in multi-modal fusion networks
• DTW-NN: A novel neural network for time series recognition using dynamic alignment between inputs and weights
• Adaptively constrained dynamic time warping for time series classification and clustering

QGIS: merge tra punti GPS e Mappa

QGIS (fino al 2013 noto come Quantum GIS) è un'applicazione desktop GIS open source che permette di visualizzare, organizzare, analizzare e rappresentare dati spaziali. Attualmente è il software GIS open source più diffuso al mondo. [fonte Wikipedia]

Nella mia esperienza mi è capitato più volte di dover elaborare mappe, dati spaziali, ... ho pian piano imparato a fare qualcosa con QGIS e vorrei condividerla con semplicissimi tutorial.

Io uso QGIS 3.4 "Madeira" (potete scaricare l'ultima versione dal sito ufficiale) e su questa release si basano gli screenshot che vedrete.

Nel post sulla dissolvenza abbiamo costruito la mappa di Milano suddivisa per ACE, ora aggiungeremo informazioni alla mappa utilizzando Openstreetmap, in particolare andremo a sovrapporre le fermate della metropolitana.

Per scaricare i punti GPS userò in particolare una libreria sviluppata da me per MagicBus, MagicWrap, un wrapper delle API Overpass scritto in R utilizzato per automatizzare alcune operazioni.
Per scaricare ed installare il package, scaricare i dati e creare un csv, eseguire i comandi di seguito (potrebbero esserci problemi di dipendenze dovuti alla versione di R, uso uguale invece di gets per problemi di formattazione del testo sul blog).

library(devtools)
install_github("MagicBusApp/magicwrap")

library(magicwrap)
milano_atm =  poi.osm(city = "Milano", tag = "railway", value = "station")

atm_to_csv = write.csv2(x = milano_atm$osm_points, "atm.csv")

Abbiamo caricato il nostro classico shapefile di Milano su QGIS.

Milano Shapefile

Ora per caricare un csv su QGIS bisogna andare sulla barra dei menù:

• Layer
• Aggiungi Layer
• Aggiungi Layer da testo delimitato

Caricamento .csv

Fermate caricate

 Ora vediamo le fermate sulla mappa, notate dei gruppi di punti, lì la fermata è stata mappata come poligono, in MagicWrap c'è una funzione che converte i poligoni in punti.
Per unire i livelli andiamo sulla barra dei menù:

• Vettore
• Strumenti di gestione dati
• Unisci gli attributi per posizione

Inserire attributi per posizione

Inserire attributi per posizione

Nella finestra selezionare come primo livello il csv e come secondo il file geografico, come operazione da eseguire l'intersezione.
Verrà creato un nuovo livello, i punti GPS che tra gli attributi riportano anche le informazioni contenute nella mappa, in questo caso quelle interessanti sono ACE e ASC.

Possiamo salvare tutto su un .csv adesso. 

Individuare pattern col SAX Encoding

Cos’è SAX Encoding?

Il Symbolic Aggregate approXimation (SAX) encoding è un metodo per semplificare le serie storiche. E’ stato inventato da Eamonn Keogh e Jessica Lin nel 2002.

SAX è un modo di trasformare una serie temporale in una sequenza di simboli. L’idea di base è che ogni simbolo rappresenta un intervallo.

Questa tecnica permette di condurre una dimensionality reduction sulle serie storiche quindi possiamo considerarlo un metodo non supervisionato.

Bisogna evidenziare che il SAX è una tecnica molto robusta ai valori mancanti e quindi utilizzabile quando si è in presenza di serie di lunghezza diversa.

Per rafforzare l’ultima considerazione vi invito a leggere un articolo di KDNuggets dove viene utilizzato SAX Encoding per trovare anomalie tra le stagioni di varie squadre di baseball. Non tutte le stagioni hanno lo stesso numero di match e alcuni dati sono mancanti.

Dataset

Visto il momento in cui viene scritto questo articolo proverò a utilizzare il SAX Encoding sui dati della protezione civile cercare i pattern tra le serie regionali dei contagi da Covid-19.

• L’orizzonte temporale và dal 24/02/2020 al 21/04/2020.
• La distribuzione geografica è per Regione anche se il Trentino - Alto Adige è diviso nelle due province autonome.

La colonna che useremo sarà totale_positivi.

Riporto direttamente dal readme del repository ufficiale della protezione civile
Totale attualmente positivi (ospedalizzati + isolamento domiciliare)

ATTENZIONE!!! IMPORTANTE.

Naturalmente non sono un medico/epidemiologo/virologo/… quindi il mio interesse è solo trovare un dataset che mi permettesse di spiegare come funziona l’algoritmo, non trarrò nessuna conclusione dai dati.

data	denominazione_regione	totale_positivi
2020-02-24T18:00:00	Abruzzo	0
2020-02-24T18:00:00	Basilicata	0
…	…	…
2020-04-21T17:00:00	Valle d’Aosta	522
2020-04-21T17:00:00	Veneto	10077

Tutto il codice di questo articolo è disponibile su un notebook Google Colab. Il dataset aggiornato può essere scaricato qui.

Come funziona l’algoritmo

Step 1: Standardizzare la serie

Per applicare questa tecnica bisogna avere la serie storica sulla riga e ad ogni colonna coinciderà uno step temporale.

Regione	1	2	…	57
Abruzzo	0	0	…	2067
Basilicata	0	…	…	245
…	…	…	…	…
Veneto	32	…	…	10077

Le serie devono essere standardizzate così da avere la stessa scala.

\(Z = \frac{x - \mu}{\sigma}\)

Questa tecnica non si basa solo sui volumi ma anche sugli andamenti. In questo caso in numeri assoluti la Lombardia è molto distante dalle altre regioni, ma lo sarà anche come andamento? Questa è la domanda a cui ci aiuterà a rispondere SAX Encoding.

tot_pos_pivot = totale_positivi.pivot(index="data", columns="denominazione_regione", values="totale_positivi").reset_index()

df_stand = ((tot_pos_pivot.drop("data", axis=1) - tot_pos_pivot.mean())/tot_pos_pivot.std()).T

Ora abbiamo la nostra matrice con sull’indice il nome della serie, ovvero la regione e sulle colonne il periodo da 0 a 57.

Regione	1	2	…	57
Abruzzo	-1.067762	-1.067762	…	1.699956
Basilicata	-1.063952	…	…	1.041431
…	…	…	…	…
Veneto	-1.244880	…	…	10766

Step 2: Piecewise Aggregate Approximation

Come accennato inzialmente questo metodo riduce la dimensionalità della serie e lo fa proprio con la P.A.A., l’idea alla base è “riassumere” la sequenza in una serie di segmenti che ne riducano la lunghezza ma ne rappresentino ugualmente l’informazione.

Data una time series \(Y =\) [\(Y_{1}\), \(Y_{2}\), …, \(Y_{n}\)] che può essere ridotta in una sequenza \(X =\) [\(X_{1}\), \(X_{2}\), …, \(X_{m}\)] con \(m \leq n\) l’equazione che sintetizza gli elementi della time series originale per adattarli alla nuova dimensione è la seguente:

\(\bar{X}_i = \frac{m}{n} \cdot \sum_{j=n/N(i-1) + 1}^{(n/M)\cdot i} x_j\)

Quindi la media per ogni finestra temporale che costituirà la nuova sequenza. Si possono notare subito due casi particolari:

• m = n, la nuova sequenza è identica alla serie storica originale
• m = 1, la nuova sequenza ha un solo valore equivalente alla media dei valori dell’intera time series

Per applicare la P.A.A. avremo bisogno di scegliere la dimensione della finestra temporale su cui calcolare l’informazione, questo influirà anche su quanti segmenti verranno utilizzati per ricostruire la serie storica.

In questo caso visto che i vari ragionamenti (sentiti in TV, ribadisco di non avere nessuna competenza) si basano su un tempo di 2 settimane imposterò il parametro (window) w=14.

Ridimensionata, ogni serie sarà una sequenza di lunghezza_serie/w periodi si prende la media del periodo come nuovo valore.

Attenzione: Bisogna impostare un controllo per non perdere informazioni, quindi se il nuovo numero di periodi non è un intero và sempre arrotondato per eccesso.

windows = 14

new_len = int(np.ceil((df_stand.shape[1]/windows)))

# Crea il nuovo dataFrame
df_PAA = pd.DataFrame(index = df_stand.index, columns = range(0, new_len))

# Calcola la media di ogni window
ind = 0
for i in range(0,df_stand.shape[1], windows):
    avg = df_stand.iloc[:,i:i+windows].mean(axis=1)
    df_PAA[ind] = avg.values
    ind +=1

denominazione_regione	0	1	2	3	4
Abbruzzo	-1.060876	-0.785615	0.352059	1.252060	1.696609
Basilicata	-1.055972	-0.846662	0.516006	1.239693	1.028541
Calabria	-1.088706	-0.797959	0.449874	1.246324	1.333274

Step 3: SAX

Preparato il dato si può lavorare a ciò che contraddistingue questa tecnica. Adesso il dato numerico verrà converito in stringa. La SAX String è derivazione del numero di livelli in cui si vuole suddividere la serie, ogni livello avrà un’etichetta, queste etichette verranno successivamente concatenate creando la SAX String. La scelta dei livelli influisce sul risultato quindi bisognerebbe sempre farsi affiancare da un esperto di dominio, il consiglio è quello di non creare troppi livelli altrimenti si avrebbe un’eterogeneità troppo alta e nono si troverebbero pattern.
I livelli possono essere fissati oppure calcolati ad ogni periodo, in questo caso si preferisce la seconda opzione.

Io non conoscendo il fenomeno supporrò 3 livelli [“A”, “B”, “C”] e suddividerò i valori:

• dal minimo fino al primo quartile
• tra primo e terzo quartile
• superiore al terzo quartile

In questi caso potremmo interpretarli come:

• pochi contagi, ottima situazione
• situazione nella media
• grave emergenza

Una volta definiti i livelli questi vengono concatenati in un’unica stringa detta appunto “SAX string”. Es. ABAAB

# Sto ridondando con i DataFrame creati ma è per mostrare i vari step, sorry.
binned = pd.DataFrame(index = df_PAA.index, columns = df_PAA.columns)

for j in range(0, df_PAA.shape[1]):
            bins = []
            bins.append(df_PAA[j].min()-.01)
            bins.append(df_PAA[j].quantile([0.25]).values[0])
            bins.append(df_PAA[j].quantile([0.75]).values[0])
            bins.append(df_PAA[j].max()+.01)
            labels = ["A", "B", "C"]
            binned[j] = pd.cut(df_PAA[j], bins, labels=labels)
binned['SAX_string'] = binned.apply(''.join, axis=1)

Il risultato finale del SAX Enconding lo si può vedere in tabella.

denominazione_regione	SAX_string
Abbruzzo	CAABC
Basilicata	CABBA
Calabria	BABBB
…	…

Come utilizzare le SAX String?

Anomaly Detection

Il SAX Encoding può essere considerato un metodo non supervisionato di anomaly detection dove per anomalia (o outlier) non si intende la singola osservazione ma l’intera serie storica.

Possiamo calcolare la frequenza di ogni SAX String, impostato il limite che per noi determina un outlier etichettare come tale tutte le serie con frequenza della SAX String minore o uguale al bound.

Nell’immagine ho impostato come bound 1

In rosso abbiamo gli outlier, blu le sequenze con frequenza maggiore di 1, vediamo che quasi ogni regione ha un suo andamento, questo dovuto a diversi fattori (territorio, controlli, …)

Clustering

Sempre rimanendo nel campo dei modelli non supervisionati le SAX String possono essere utilizzate per il clustering.

Un’applicazione che mi viene in mente è il clustering gerarchico utilizzando come misura di similarità (o dissimilarità a seconda se agglomerativo o divisivo) la distanza di Levenshtein.

Classification

Passando ai metodi supervisionati, avendo già delle serie storiche etichettate, penso ad esempio a quelle sul consumo energetico che hanno anche una certa stagionalità, si potrebbe applicare il KNN alla sequenza per classificare le nuove serie.

Stream Data

Nel precedente paragrafo ho citato la stagionalità perchè in questo articolo non è stato trattato un tema di forte interesse come il dato in stream.
Con la diffusione dell’IoT le time series in stream sono un tema sempre più caldo sia per quanto riguarda l’immagazzinamento dei dati che per le analisi.

Nel caso di serie con stagionalità, o forti componenti cicliche, oltre alla finestra temporale per la P.A.A. si potrebbe pensare ad una finestra per l’aggiornamento dei periodi con lo stesso timing dell’effetto. Ad esempio, se ho una serie con ciclicità mensile e voglio valutare i pattern su finestre settimanali potrei impostare:

• orizzonte di aggiornamento mensile
• finestra per il SAX Encoding settimanale

così da studiare SAX String di lunghezza 4 o 5.

Naturalmente con fenomeni fortemente erratici questo è più complicato da decidere, quasi sicuramente i risultati dei nostri modelli (anomaly detection, classificazione, clustering, …) varierebbero molto frequentemente.

Prendendo ad esempio il nostro dataset. Se ci fermassimo una settimana prima, nel grafico seguente i dati si fermano al 13 aprile, le sequenze non singleton sarebbero molte di più. Questo solo per rimarcare ancora una volta l’importanza della scelta dela giusta finestra temporale e del periodo di osservazione in relazione al fenomeno che si và ad analizzare.

Se volete approfondire questo tema vi consiglio il paper “A Symbolic Representation of Time Series, with Implications for
Streaming Algorithms” di Jessica Lin, Eamonn Keogh, Stefano Lonardi, Bill Chiu.

QGIS: calcolare un'area

QGIS (fino al 2013 noto come Quantum GIS) è un'applicazione desktop GIS open source che permette di visualizzare, organizzare, analizzare e rappresentare dati spaziali. Attualmente è il software GIS open source più diffuso al mondo. [fonte Wikipedia]

Nella mia esperienza mi è capitato più volte di dover elaborare mappe, dati spaziali, ... ho pian piano imparato a fare qualcosa con QGIS e vorrei condividerla con semplicissimi tutorial.

Io uso QGIS 3.4 "Madeira" (potete scaricare l'ultima versione dal sito ufficiale) e su questa release si basano gli screenshot che vedrete.

Partendo dal layer costruito nel tutorial precedente aggiungiamo una proprietà ai nostri dati calcolando l'area di ogni ACE.

Per prima cosa bisogna andare nella tabella attributi e renderla modificabile dall'icona a forma di matita in alto a sinistra.

Modifica Tabella Attributi

Abilitare la creazione di un "nuovo campo calcolabile".

Nuovo Campo Calcolabile

Dalla maschera impostiamo il nome del campo in uscita e il tipo, in questo caso:

• nome: AREA_KMQ
• tipo: numerico (real) con 3 decimali

Dal menù delle funzioni scegliamo Geometria --> $area.

• $area misura l'area dalle coordinate geografiche (non si vedono nella tabella attributi)
• area misura l'area degli attributi
• l'area è misurata in metri quadrati quindi poi dovremo convertirla ( $area/1000000)

Calcolo dell'area delle ACE

 In coda alla tabella è stato aggiunto il nostro campo calcolato.

Campo calcolato

Fatto questo possiamo salvare ed uscire dalla modalità di modifica della tabella attributi.

QGIS: dissolvenza

QGIS (fino al 2013 noto come Quantum GIS) è un'applicazione desktop GIS open source che permette di visualizzare, organizzare, analizzare e rappresentare dati spaziali. Attualmente è il software GIS open source più diffuso al mondo. [fonte Wikipedia]

Nella mia esperienza mi è capitato più volte di dover elaborare mappe, dati spaziali, ... ho pian piano imparato a fare qualcosa con QGIS e vorrei condividerla con semplicissimi tutorial.

Io uso QGIS 3.4 "Madeira" (potete scaricare l'ultima versione dal sito ufficiale) e su questa release si basano gli screenshot che vedrete.

Ripartendo dal post precedente abbiamo i dati del comune di Milano.

Comune di Milano

La mappa ha tutti i livelli che come si può vedere dalla tabella attributi sono parecchi.

Attributi Shapefile

Se noi volessimo semplificare la mappa mantenendo ad esempio solo il livello ACE visibile ci basterebbe utilizzare uno degli Strumenti di Geoprocessing messi a disposizione da QGIS, la DISSOLVENZA.

Dalla barra dei menù selezionare:

• Vettore
• Stumenti di Geoprocessing
• Dissolvi

Selezionare per quale attributo si vuole dissolvere la mappa, in questo caso le ACE.

Selezione Attributo

 Eseguita la dissolvenza verrà creato un layer che mostra come livello principale gli attributi per cui abbiamo dissolto, tutti gli altri sono stati inglobati a livello visivo ma restano nella Tabella Attributi.

Mappa delle ACE di Milano

QGIS: filtrare un'area

QGIS (fino al 2013 noto come Quantum GIS) è un'applicazione desktop GIS open source che permette di visualizzare, organizzare, analizzare e rappresentare dati spaziali. Attualmente è il software GIS open source più diffuso al mondo. [fonte Wikipedia]

Nella mia esperienza mi è capitato più volte di dover elaborare mappe, dati spaziali, ... ho pian piano imparato a fare qualcosa con QGIS e vorrei condividerla con semplicissimi tutorial.

Io uso QGIS 3.4 "Madeira" (potete scaricare l'ultima versione dal sito ufficiale) e su questa release si basano gli screenshot che vedrete.

Per questo esempio verrà utilizzato uno shapefile della Lombardia scaricato dal sito dell'ISTAT.

Schermata Avvio QGIS

Per caricare il layer basta cercarlo tra le nostre cartelle, è un file .zip, entrare e trascinare nel riquadro in basso il file .shp. Una volta caricato ci apparirà la mappa nel riquadro centrale.

Shapefile Caricato

Ora il nostro obiettivo è isolare un territorio, ad esempio il comune di Milano.

Per prima cosa facciamo click col tasto destro sul layer caricato e scegliamo "Apri tabella attributi".

Apri tabella attributi

Tabella attributi

Per filtrare bisogna scegliere l'icona (ho provata a cerchiarla in rosso nell'immagine precedente) con una piccola epsilon, faccio notare che il file ISTAT non è modificabile.
Se diamo un'occhiata alle colonne notiamo che ci sono un sacco di informazioni/livelli in questo file. Infatti solo rimanendo sui confini amministrativi abbiamo:

1. Confini regionali
2. Confini provinciali
3. Confini comunali
4. Confini subcomunali (ASC - area sub comunale)
5. Confini censimento (ACE - area censimento)

Indichiamo che vogliamo filtrare per codice ISTAT, volendo il comune, il codice ISTAT di Milano è 003015146:

• 003 Lombardia
• 015 Provincia di Milano
• 146 Comune di Milano

Nel caso del file il codice è un numerico quindi "cadono" i primi 2 zero. Il codice per la selezione è "COD_ISTAT" = 3015146

Selezione comune

Click su "Seleziona Elementi" e avremo filtrato l'area di interesse che cambierà colore sulla mappa.

Selezione Milano

Trattare i valori mancanti nelle serie storiche

Introduzione

Il trattamento e l’imputazione dei valori mancanti (missing values) è uno step molto delicato per ogni progetto di data science.

Esistono diverse strategie per l’imputazione e tutte possono portare a errori perchè si sta introducendo un dato “artificiale”.

Un consiglio che viene dato spesso è, in fase di imputazione di valori mancanti creare per ogni feature che si tratta una nuova variabile booleana “nomeFeature_isMissing” per tracciare quali valori sono reali e quali indotti durante il processo di cleaning.

Solitamente gli step che vengono seguiti in questa fase del preprocessing del dataset sono i seguenti:

• Se la percentuale di valori mancanti è alta (la soglia varia a seconda del contesto) non si può considerare la variabile, quindi viene eliminata uan feature
• Non eliminare mai l’intera osservazione a meno che non abbia valori mancanti su ogni feature
• Scegliete una tecnica di imputazione basandovi sul tipo di dato e fenomeno che si sta trattando

Lo strumento Python più utilizzato per questo compito è il SimpleImputer di scikit-learn che offre quattro strategie di sostituzione:

1. sostituire i missing values con la media
2. sostituire i missing values con la mediana
3. sostituire i missing values con la moda
4. sostituire i missing values con una costante

C’è una cosa da osservare, nessuna delle strategie descritte si adatta bene allo studio delle serie storiche. Per dimostrare questa affermazione faremo un test utilizzando la serie storica Airpassenger, una sequenza mensile del numero di passeggeri sui voli internazionali tra il 1949 e il 1960.

Tutto il codice di questo articolo è disponibile su un notebook Google Colab.

Creazione dei Missing Values

Airpassenger non presenta valori mancanti quindi verranno introdotti in modo completamente casuale.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

airpassengers = pd.read_csv("airpassenger.csv", index_col = 0 )
airpassengers.columns = ["Passengers"]
airpassengers.head()

Month	Passengers
1949-01	112
1949-02	118
1949-03	132
1949-04	129
1949-05	121

Verrano generati 10 indici mediante la funzione randint del modulo random di numpy. In questi dieci indici verrà sostituito il valore originale con NaN.

NOTA. Visto come vengono generati i valori mancanti il risultato cambierà ogni volta che viene eseguito il codice

airpassengers_MV = airpassengers.copy()
na_index = np.random.randint(0, airpassengers.shape[0]-1, 10)
airpassengers_MV.iloc[na_index,0] = np.nan

plt.figure(figsize=(16,8))
plt.plot(range(airpassengers.shape[0]), airpassengers_MV["Passengers"].values, color='tab:blue')
plt.gca().set(title="Airpassengers - Missing Values", xlabel="Time", ylabel="Passengers")
plt.show()

Month	Passengers
1949-01	112
1949-02	118
1949-03	132
1949-04	NaN
1949-05	121

Imputazione dei Missing Values

Di seguito verranno testate sulla serie tutte le strategie proposte da SimpleImputer di Scikit-learn.

from sklearn.impute import SimpleImputer

SimpleImputer rispetta la classica interfaccia di sklearn, una volta inizializzato l’oggetto può essere utilizzando sfruttando i metodi di base, quali fit, transform. Nel nostro esempio non avendo interesse verso il forecasting e quindi non dividendo la serie in train e test, useremo fit_transform su tutta la sequenza.

Il metodo non è altro che la concatenazione dei due:

• fit calcola il valore da imputare (media, mediana, più frequente) sulla serie
• transform applica realmente l’imputazione

In input necessita di un oggetto a 2D, noi passeremo un DataFrame, ma restituisce un array numpy.

Sostituzione con la media

La strategia di SimpleImputer viene applicata mediante il parametro strategy, come detto in precedenza si può scegliere tra:

• “mean” default
• “median”
• “most_frequent”
• “costant”

imputer_mean = SimpleImputer(strategy = "mean")
airpassengers_mean = imputer_mean.fit_transform(airpassengers_MV)

airpassengers_mean[:5]

array([[112. ],
[118. ],
[132. ],
[282.23880597],
[121. ]])

in grassetto il valore imputato

Questo approccio risente molto del trend, Airpassenger è una serie con trend lineare crescente e un effetto stagionale moltiplicativo (le onde sono sempre più alte ogni anno), l’approccio dell’imputazione utilizzando la media crea picchi all’inizio della sequenza e valli sul finire, anche se queste ultime sono molto meno evidenti dei primi. Se la serie fosse ancora più lunga questo fenomeno sarebbe sempre più marcato. Di seguito il confronto con l’originale (in rosso).

Sostituzione con la mediana

La mediana è una misura di sintesi robusta agli outlier, quindi si potrebbe pensare che l’effetto trend venga un pò mitigato utilizzando questa strategia.

airpassengers_median = imputer_median.fit_transform(airpassengers_MV)
airpassengers_median[:5]

array([[112.],
[118.],
[132.],
[268.],
[121.]])

L’unica cosa che è cambiata è che la mediana è più bassa della media, quindi i picchi sono più bassi ma si acuisce l’effetto nella parte finale della serie.

Le strategie con media e mediana possono essere chiaramente scartate per l’imputazione di missing values nelle serie storiche

I paragrafi successivi servono solo per testare anche le restanti opzioni presenti tra le strategie di SimpleImputer. Solitamente queste vengono utilizzate su variabili categoriche non numeriche continue. Quindi sappiamo a priori che non si adatteranno bene alle Time Series.

Sostituzione con la moda

imputer_mode = SimpleImputer(strategy="most_frequent")
airpassengers_mode = imputer_mode.fit_transform(airpassengers_MV)
airpassengers_mode[:5]

array([[112.],
[118.],
[132.],
[229.],
[121.]])

Visto che il trend cresce linearmente e la stagionalità con effetto moltiplicativo la moda è un valore più basso della media, ma anche della mediana. In questo caso sembra dare un andamento quasi naturale alla serie nella parte centrale, ma quando il valore mancante si trova nei primi/ultimi periodi l’errore sistematico è evidente.

N.B. La strategia, “Most Frequent”, è molto utile se si hanno missing in variabili categoriali come detto in precedenza. Un’altra strategia in questi casi può essere applicare un modello di Machine Learning (es. KNN) usando come target la variabile in cui dobbiamo imputare i missing values e come features tutte le altre, creando il training set per allenare il modello sui non missing e applicarlo sul test, ovvero le osservazioni che presentano i dati mancanti

Sostituzione con una costante

Il primo problema che ci si pone in questo caso è “Quale costante utilizzare?”:

• un quantile?
• il minimo/massimo?
• un numero a caso?

Per questo test verranno usati lo zero, il minimo e il massimo.

Sostituzione con zero

Lo zero solitamente è sconsigliatissimo perchè non si capisce se è assenza di valore o presenza pari a zero, crea molta abiguità.

Esempio: presenza ad eventi, se viene effettuata la sotituzione con zero non si riesce più a riconoscere se è un dato imputato o effettivamente c’è presenza zero (in aggiunta al restante errore indotto).

imputer_zero = SimpleImputer(strategy="constant", fill_value=0)
airpassengers_zero = imputer_zero.fit_transform(airpassengers_MV)

Sostituzione con massimo/minimo

Queste due strategie ci daranno Time Series completamente sballate (dedicheremo pochissime righe)

imputer_max = SimpleImputer(strategy="constant", fill_value=airpassengers_MV.Passengers.max())
airpassengers_max = imputer_max.fit_transform(airpassengers_MV)

Analogamente per il minimo

imputer_min = SimpleImputer(strategy="constant", fill_value=airpassengers_MV.Passengers.min())
airpassengers_min = imputer_min.fit_transform(airpassengers_MV)

Imputazione con Media Mobile

La strategia vincente con le serie storiche è quella dell’imputazione con la media mobile.
Esistono diversi tipi di media mobile (qui la [pagina Wikipedia] (https://it.wikipedia.org/wiki/Media_mobile)), nella maggior parte dei casi può bastare quella lineare che verrà utilizzata in questo articolo.

Data un serie storica ${y_t}$ con $t=(1, 2, ..., T)$, sia un generico elemento della serie $t$ un valore mancante e data una finestra temporale di dimensione $N$, siano:

• $m_1$ gli $N$ periodi antecedenti il valore mancante
• $m_2$ gli $N$ periodi successivi al valore mancante
• $\theta_i$ il peso da attribuire all’i-esimo valore osservato. (Per noi sarà pari ad 1 visto che vogliamo una media aritmetica semplice).

Si definisce media mobile al tempo $t$:

$mm_i$ = $\frac{1}{k}$ $\sum_{i=-m_1}^{m_2}$ $\theta_i$ $y_{t+1}$

Dove $k = {m_1} + {m_2} + 1$

NOTA. In questo articolo non useremo la funzione rolling_mean di Pandas per mostrare come funziona la media mobile ma verrà costruita un piccolo script ad-hoc.

L’imputazione con questa tecnica permette di imputare al valore mancante la media locale in un range deciso dall’analista a seconda del fenomeno che si và a studiare.

Il metodo è un pò più esoso in termini di calcolo e presenta alcuni problemi che vanno risolti, ad esempio:

• se manca il valore all’inizio o alla fine della serie?
• se nel range indicato ci sono più valori mancanti?
• se abbiamo missing values contigui?

Sta al data scientist scegliere quale soluzione applicare ad ogni domanda a seconda del contesto.

Nel nostro caso applicheremo una finestra temporale tra i tre valori prima e i tre valori dopo il missing value, gestendo i due casi:

1. Il lower bound è negativo
2. L’upper bound supera la lunghezza della serie.

Se nel range si presentano più missing values viene semplicemente calcolata la media tra i valori presenti.

Esempio:
[2, NaN, 3] --> [2, 2.5, 3]
[2, NaN, NaN] --> [2, 2, 2]

La prima cosa che si nota da questo esempio è che cade il +1 finale nel calcolo di k visto che quel valore per noi è mancante.

airpassengers_MA = airpassengers_MV.copy()
steps = 3

for idx in na_index:
    lower = idx - steps
    upper = idx + steps + 1
    if lower<0:
        lower=0
    if upper>airpassengers_MA.shape[0]:
        upper=airpassengers_MA.shape[0]
    
    airpassengers_MA.iloc[idx,0] = airpassengers_MA.iloc[lower:upper,0].mean()

L’andamento sembra abbastanza naturale, questo è confermato dal confronto con la serie storica originale. Quindi per le serie storiche la media mobile è nettamente la strategia vincente.

Imputare un valore mancante introduce sempre errore sistematico, si può solo scegliere la strategia più adatta per minimizzarlo.