Selezioni ibride su un DataFrame pandas multi-indice

La selezioni di viste da un dataframe è cosa basilare per un Data Scientist. L’altro giorno al lavoro ho però incontrato una situazione spinosa:

selezionare dati da un dataframe multi-indice con maschera booleana e con specifico valore in indice.


La soluzione è il metodo pandas.DataFrame.query()

Vediamo un esempio pratico. Il nostro input è:

import pandas as pd
index = pd.MultiIndex.from_product([['xy1','xy2','xy3'], ['1','2','3','4','5']], names=['Tag', 'Page'])
df = pd.DataFrame([1,1,1,4,5,1,1,61,4,51,1,1,4,5,1], index, columns=['Value'])

mentre ciò che vogliamo ottenere è questo:

index2 = pd.MultiIndex.from_product([['1','2','3']], names=['Page'])
df2 = pd.DataFrame([1,1,1], index2, columns=['Value'])
df2

Usando pandas.DataFrame.query() possiamo risolvere il dilemma in molto molto pythoniano:

df1 = df.query("Tag == 'xy1' & Value == 1").reset_index(level=0, drop=True)

(il .reset_index() è utilizzato per resettare l’indice e ottenere il risultato cosmetico voluto)

 

Forecasting di vendite: serie storiche e ML

Predirre le vendite è il sogno più bagnato di ogni manager che si rispetti. Difatti, avere stime affidabili sull’andamento del mercato è la chiave del budgeting: più le stime sono corrette, più la programmazione diventa facile ed effettiva.

Come in ogni problema di predizione, la variabile che influisce di più sulla bontà del modello è l’informazione disponibile. Nella mia esperienza lavorativa, la quasi totalità delle aziende ha (per ora) una cultura del dato molto scarsa, se non nulla. Il data scientist si troverà dunque a dover lavorare con dati parziali, da reclamare tramite molteplici fonti, la cui affidabilità sarà semplicemente dubbia. Ovviamente, tale ostacolo è prettamente individuale e varierà da scenario a scenario: è dunque fuori dallo scopo di questo post.

Quello di cui voglio veramente parlare è come approcciare the algoritmic side of thigs . Ordunque, vi presento qui uno scoppiettante programmino per affrontare il forecasting delle serie storiche come dei veri pro.


Le basi

Il problema può essere affrontato in base a due variabili interdipendenti:

  1. La finestra di predizione
  2. Il lag delle variabili

La finestra di predizione rappresenta con quanto anticipo vogliamo predirre Y (e.g. se la predizione ha 3 mesi di anticipo, la finestra sarà di 3 mesi). Il lag delle variabili rappresenta invece la relazione temporale con cui le variabili indipendenti influenzano la dipendente. Ad esempio: le vendite di un’industria metallurgica saranno influenzate dall’andamento del PIL globale. Ma questa relazione non sarà immediata temporalmente: l’aumento di domanda del metallo sarà successivo di x mesi all’aumento del PIL.

Mentre la finestra di predizione è decisa dal Data Scientist e dal management, il lag delle variabile è fondamentalmente ignoto. Per trovarlo, in maniera ottimale, si ricorre alla Cross Validation. L’approccio che io propongo è univariato.

Restate con me, così vi spiego le funzioni per questo scopo.

L’oracolo algoritmico

Il programma è costituito da tre funzioni:

  1. laggedDataMat: lagga il dataset (ossia “sposta” le variabili di un determinato lag)
  2. bestLagRegr: trova il lag ottimale (univariato) per ogni variabile
  3. modelSelection: prende l’outup di bestLagRegr ed esegue una grid search per trovare l’ottimale modello predittivo

Le librerie richieste per l’esecuzione sono:

import numpy as np, pandas as pd, pickle, copy
from sklearn.model_selection import cross_val_score, TimeSeriesSplit, GridSearchCV
from sklearn import ensemble
from sklearn.ensemble import GradientBoostingRegressor

Ecco qui di seguito le funzioni in Python:

laggedDataMat
def bestLagRegr(dataMat, maxLag, minLag, yName):

    # Purpose: identify the best lag for variables to be used in regression. Uses CrossValidation
    # (folds TimesSeriesSplit)

    # data - a panda DF with named columns
    # maxLag - the maximum lag to be tested. Integer
    # minLag - the minimum lag to be tested. Integer
    # yName - the name of the column of 'data' containing the dependent var. String

    data = dataMat.copy()
    colnames = [y for y in data.columns if y != yName]
    lags = range(minLag, maxLag)
    folds = TimeSeriesSplit(n_splits=int(np.round(data.shape[0] / 10, 0)))

    results = {}

    for col in colnames:

        scores = []
        lags_list = []

        for l in lags:
            varname = col + str('_lag') + str(l)
            data[varname] = data[col].shift(l)
            YX = data[[yName, varname]]
            YX = YX.dropna().as_matrix()

            # Build regressor and estimate metric of prediction performance with CV
            regr = ensemble.GradientBoostingRegressor(learning_rate=0.01, max_depth=1, n_estimators=500)
            perform = cross_val_score(regr, X=YX[:, 1].reshape(-1, 1), y=YX[:, 0], cv=folds,
                                      scoring='neg_median_absolute_error')

            # Store scores result and lags
            scores.append(np.median(perform))
            lags_list.append(l)

        # Calculate best score, corresponding best lag and store it in a dictionary, containing all colnames
        best_score = max(scores)
        best_lag = lags_list[scores.index(best_score)]
        results[col] = [best_lag, best_score]

    return(results)
bestLagRegr
def laggedDataMat(dataMat, yName, lagsDict):

    # Purpose: build a lagged DF
    # dataMat: the unlagged DF with named columns
    # yName: name of the dependent var. String
    # lagsDict: dictionary produced with 'bestLagRegr, containing:
    #                                                               - keys: column names of dataMat
    #                                                               - elements: lists with lag, CV score
    #
    # Output: a panda DataFrame with columns order sorted alphabetically

    # Initialize empty DF
    df = pd.DataFrame(index=dataMat.index)
    # Set dependent var
    df[yName] = dataMat[[yName]]

    # Creating and adding the lagged vars
    for colName in lagsDict.keys():
        l = lagsDict[colName][0]
        colNameLag = colName + str('_lag') + str(l)
        df[colNameLag] = dataMat[[colName]].shift(l)

    df = df.sort_index(axis=1)

    return(df)
modelSelection
def modelSelection(maxLag, data, depVar, toSave, pSpace = None, alpha = 0.95):

    if pSpace is None:
        pSpace = dict(n_estimators=list(range(5, 2000, 10)),
                      learning_rate=list(np.arange(0.001, 1, 0.1)),
                      max_depth=list(range(1, 3, 1)))

    lags = range(1, maxLag)
    results = dict()

    for lagMin in lags:
        print('Esimating model for lag: ', lagMin)
        lagAnalysis = bestLagRegr(data, maxLag, lagMin, depVar)
        lagMat = laggedDataMat(data, depVar, lagAnalysis)
        lagMat = lagMat.dropna()
        trainY = np.ravel(lagMat[depVar])
        lagMat = lagMat.drop([depVar], 1)

        folds = TimeSeriesSplit(n_splits=int(round(lagMat.shape[0] / 10, 0)))

        model = GradientBoostingRegressor(loss='ls')
        regr = GridSearchCV(estimator=model, param_grid=pSpace, scoring='neg_mean_squared_error', cv=folds)
        regr.fit(lagMat, trainY)

        modelName = toSave + '/' + 'modelOI' + '_lag' + str(lagMin) + '.sav'
        pickle.dump(regr.best_estimator_, open(modelName, 'wb'))
        temp = dict(BestModel=regr.best_estimator_, Score=regr.best_score_, Lags=lagAnalysis)

        regrQUpper = copy.deepcopy(regr.best_estimator_)
        regrQUpper.set_params(loss='quantile', alpha=alpha)
        regrQUpper.fit(lagMat, trainY)
        temp['QUpper'] = regrQUpper

        regrQLower = copy.deepcopy(regr.best_estimator_)
        regrQLower.set_params(loss='quantile', alpha=(1-alpha))
        regrQLower.fit(lagMat, trainY)
        temp['QLower'] = regrQLower

        key = 'Lag ' + str(lagMin)
        results[key] = temp

    return(results)

Le Macchine possono vederti

Sulla costruzione di algoritmo di riconoscimento visivo


Avrete certamente notate la mia più che settimanale assenza. Purtroppo sono stato impegnato per motivi lavorativi: per farmi perdonare, ho deciso di produrre un bell’articolo (con tanto di programma Python) su un argomento iper attuale: la Machine Vision.

Intro

La Machine Vision raccoglie tutte le tecniche algoritmiche per il riconoscimento automatico di immagini. In essenza, utilizzando queste tecniche, è possibile insegnare ad una macchina a riconoscere gli oggetti più disparati (anche in tempo reale): l’applicazione più semplice é ad esempio la lettura automatica di testi. Come quando vi arriva la multa a casa, perché siete entrati come dei fessi nella ZTL: la telecamera, partendo da una foto, ha automaticamente letto la targa del vostro carro.

La teoria

A prima vista potrebbe sembrare che algoritmi di tale portata possano essere complessi e sostanzialmente distaccati dalla teoria di Data Science. In realtà, l’algebra lineare è ancora una volta onnipresente e viene in aiuto. Difatti, i metodi di Machine Vision non differiscono di molto dai metodi di Data Science in generale. Le differenze più marcate sono in come i dati (i.e. le immagini) vengono pre-processate, prima di essere inviate ad un classificatore (e.g. una Random Forest o una Rete Neurale).

Ordunque, qual’è la logica generale da seguire? Le immagini presentano un inevitabile problema di fondo: sono “pesanti” dal punto di vista di utilizzo memoria. Infatti, anche una piccola foto in b\n di dimensioni (in pixel) 100×100 rappresenta, per un calcolatore, una matrice di 10.000 elementi. Triplicate il valore per una in RGB.

Dunque, a meno di avere a disposizione una elevata potenza di calcolo, il classificatore impiegherà un tempo non trascurabile per stimare tutti i parametri.

Per ovviare a questo problema, non si deve nutrire il modello con le immagini pure, ma bisognerà estrarre delle misure (in inglese features) che riassumano l’informazione contenuta in esse. Di queste misure ne esistono di tutti i gusti ed è sostanzialmente impossibile riassumerle tutte. Come esempio, si prendano quelle utilizzate nel programma descritto più avanti nell’articolo.

Le features diventano dunque le nostre nuove variabili. In linea generale, il nostro problema viene nuovamente ridotto al generale:

Y=f(X)

dove X è un vettore di lunghezza p contenente i valori per le p misure e  Y= \in \{0,1, . . ., k\} rappresentante le k classi di immagini.

Un semplice esempio per chiarire: costruiamo un algoritmo per permettere al computer di riconoscere il colore blu dal rosso e viceversa, utilizzando immagini di 10×10 pixels. In questo caso avremo n immagini, alcune completamente rosse e altre blu. Per abbiamo scelto 120.

Innanzitutto, importiamo le librerie che ci serviranno:

from PIL import Image
import numpy as np
import os
from os import walk
from sklearn.ensemble import RandomForestClassifier as RF
from sklearn import cross_validation

Ora doppiamo si procede con la creazione delle cartelle e delle immagini. Al tempo stesso si deve anche creare il vettore della variabile dipendente (indicato con Y).


# -------------------------------------------------------------------------------------- #
# Controlla l'esistenza di due cartelle per salvare le immagini. Se non presenti, le crea
directories = ["/foo/rosso/", "/foo/blu"]
for dir in directories:
    if not os.path.exists(dir):
        os.makedirs(dir)
# -------------------------------------------------------------------------------------- #

files = []
t = 0  # contatore
size = 10 # dimensione k, in pixel, di immagini kxk
iterations = range(0, 60) # costruisce 60 immagini
Y = []

for u in iterations:

    k1 = int(np.random.uniform(0, 100)) # Valori aleatori per introdurre variabilità nei colori
    k2 = int(np.random.uniform(0, 100))
    k3 = int(np.random.uniform(0, 100))
    k4 = int(np.random.uniform(0, 100))

    col = [["rosso/", (255, k1, k2)], ["blu/", (k3, k4, 255)]] # paletta colori

    for n,i in col:
        im = Image.new("RGB", (size, size))
        pix = im.load()

        for x in range(size):
            for y in range(size):
                pix[x, y] = i # setta il colore dei pixel
                
        filename = directory + n + str(t) + ".png"
        im.save(filename, "PNG")
        files.append(filename) # crea la lista dei filenames
        Y.append(n[0:-1]) # Crea il vettore variabile dipendente con le classi
        t += 1

Una volta creato il dataset e la variabile dipendente, si può procedere all’estrazione delle misure, che andranno a creare la nostra matrice X delle variabili indipendenti. Come già accennato, la scelta di quali features includere è notevolmente specifica al problema da risolvere. In questo semplice esempio, la distribuzione dei colori dei pixel sui 3 canali (RGB) è sufficiente.

Nota tecnica: nell’ultimo passaggio logico, se in R, sarebbe meglio vettorizzare la funzione anziché evocarla in loop. Ma usando Py, il for non crea troppi problemi.

print("Extracting features")
def extract_features(path):

    ''' Extract features from image of MetroBank Clapham Junction: colours distributions, sdev of the colours, number
    of objects and several measurements on those objects.
    Input: path - string with path of the picture
    Output: a vector '''


    # Resizing to 1000x1000, pixel distribution for each colour channel according to a fixed number of bins
    im = Image.open(path)
    im_arr = np.array(im)

    hist_R = np.histogram(im_arr[:, :, 0].flatten()) # Distibuzione Rosso
    hist_G = np.histogram(im_arr[:, :, 1].flatten()) # Distribuzione Verde
    hist_B = np.histogram(im_arr[:, :, 2].flatten()) # Distribuzione Blu

    features = np.concatenate((hist_R[0], hist_G[0], hist_B[0], hist_R[1], hist_G[1], hist_B[1]), axis=0)

    features = np.reshape(features, (1, len(features)))
    return (features)

# ------------------------------------------------------------------------------- #
# Evoca la funzione 'extract_features()' in loop. In R sarebbe meglio vettorizzare.

X = np.zeros((len(files), 63)) # ncol da modificare a seconda del'entrata
i = 0
for im in files:
    feat = extract_features(im)
    X[i, 0:len(feat[0])] = feat[0]
    i += 1
    print(i/len(files)*100, '%', ' done')

E alla fine, basta costruire un classificatore: in questo caso, data la semplicità dell’esempio, ho scelto una normale Random Forest. Facile da utilizzare avendo pochi parametri, non crea troppi problemi di overfitting e funziona straight out of the box.

# -------------------------------------------------------------------------------------------- #

# Stima del modello: Random Forest

print("Training")
# n_estimators is the number of decision trees
# max_features also known as m_try is set to the default value of the square root of the number of features
clf = RF(n_estimators=100, n_jobs=3)
scores = cross_validation.cross_val_score(clf, X, Y, cv=5, n_jobs=1) # Effettua la Cross Validation con 5 folds
print("Accuracy of all classes")
print(np.mean(scores)) # Stampa la media del punteggio di Cross Validation (accuracy in questo caso)

Se eseguite il programma, vedrete che l’accuracy media è di 1: abbiamo dunque costruito un modello perfetto per distinguere i colori di semplici immagini. L’accuracy di 1 tuttavia non è sempre un buon segno: potrebbe voler dire che il nostro modello è in over-fitting, ossia è troppo adattato ai dati osservati ma non è generalizzabile. In questo caso però non rischiamo di cadere in questo problema poiché il training set é molto semplice e praticamente privo di rumore.

Enjoy!