Blog

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)

Calcolo stock minimi con R

Stimare lo stock minimo delle materie prime è una procedura difficile e suscettibile ad errori, soprattutto se eseguita con metodi euristici.

Per fortuna, la statistica ci salva ancora una volta, assicurandoci un livello ottimale di stock e massimizando i risparmi dell’inventario.

Il problema

Come si presenta il problema? Abbiamo 3 variabili da considerare:

  • Tempi di riordino (in questo esempio, 10 giorni di lead time)
  • Il consumo storico
  • Il Lean Manager che ti rincorre nel panico

Per l’esempio di questo post, il prodotto da stoccare saranno fogli di lamiera (misuarti in Kilogrammi), con un tempo di riordino di 10 giorni.

Exploratory Data analysis
Fig. 1 Grafico consumo lamiera per la variabile temporale. La stagionalità è marcata

Un semplice grafico (Fig. 1) del consumo lamiera contro le decine dei giorni ci presenta una stagionalità abbastanza marcata: sarebbe dunque opportuno ottimizzare seguendo la variabile temporale, in modo da risparmiare sull’inventario.

Fig. 2 Istogramma del consumo lamiera. L’approssimazione a distribuzione Gamma è chiaramente visibile.

Se eseguiamo un istogramma del consumo lamiera, la distribuzione è inoltre chiaramente approsimabile da una Gamma.
Siamo fortunati: possiamo già azzardare un modello.

Il modello

Come creare una stima per gli stock minimi? Un problema simile è facilmente risolvibile con il concetto di intervallo di confidenza (IC).

Possiamo quindi calcolare il consumo medio di lamiera data la variabile temporale, ossia la soluzione al problema risulta costruire un intervallo di confidenza (IC) intorno a:

E[Y|(X,\omega)] = \mu_{Y|(X, \omega)}

Dove:

  • è il consumo di metallo
  • X è la matrice delle variabili indipendenti (i.e. l’indice del periodo di 10 giorni)
  • ω è il codice prodotto (in questo caso solo un tipo di lamiera)
  • μ indica la media

Per stimare questa media, è necessario fittare un GLM con famiglia Gamma e link-function logaritmica:

\Hat{y}_{i,k} = e^{\beta_{0} + \beta_{1}x_{i} + \beta_{2}\omega_{k}}

Ricordo al caro lettore che a livello della link-function, la realtà è Normale. Dato ciò, costruire un IC α è banale: bisogna solamente eseguire il classico calcolo con la distribuzione t (perchè la dev std della popolazione non è nota, mascalzone!)

Il GLM (e il così costruito IC), ci stimano il consumo medio per periodo: dunque, il limite massimo di tale IC risulta essere la scelta migliore per lo stock minimo.

Ovviamente, la decisione del Lean Manager risulta ora concentrata sulla scelta del parametro α.

Ed ecco a voi il risultato:

Stock minimi lamiera: il livello minimo suggerito è rappresentato dalla linea rossa (95% di confidenza), mentre i punti verdi sono il consumo medio per periodo.

Per una trattazione più matematicamente rigorosa del problema, potete riferirvi al mio GitHub, dove troverete la documentazione del progetto. Qui in calce, invece, troverete il codice usato.

set.seed(123)


# Prototyping for Minimum stock of metal sheets

lamiera <- read.csv2("foo/Lamiera.csv", sep= ";", dec = ",")
lamiera <- lamiera[lamiera$Year >= '2014',] # Only considering years 2014 onwards. By EDA boxplot(quantity~Year), 2014 is the first year of greater variability

# Fitting model: glm, Gamma distr with log-link function
lmTest <- glm(Quantity~as.factor(TenDaysOfYear)+as.factor(PartNo), data = lamiera, family= Gamma(link= 'log'))


# Estimating CI and plotting

alpha <- 0.98 # Significance level

codiciLamiera <- unique(lamiera$PartNo)
for(codLamiera in codiciLamiera){
  
  
  preddata <- data.frame('TenDaysOfYear' = seq(1,36, by=1), 'PartNo' = codLamiera)
  preds <- predict(lmTest, newdata = preddata, type = "link", se.fit = TRUE)
  
  critval <- qt(alpha, nrow(lamiera)-1)
  upr <- preds$fit + (critval * preds$se.fit)
  lwr <- preds$fit - (critval * preds$se.fit)
  fit <- preds$fit
  fit2 <- lmTest$family$linkinv(fit)
  upr2 <- lmTest$family$linkinv(upr)
  lwr2 <- lmTest$family$linkinv(lwr)
  
  name <- paste('foo/Lamiera.csv','lamiera_', codLamiera,'.png', sep = '')
  
  png(filename = name)
  plot(Quantity~TenDaysOfYear, data= lamiera[lamiera$PartNo == codLamiera,], main=codLamiera)
  lines(upr2, col="red")
  points(fit2, col="green")
  dev.off()
  
  name <- paste('foo/Lamiera.csv','lamiera_', codLamiera,'.csv', sep = '')
  preddata$PredictedConsumption <- round(fit2, 2)
  preddata$MinimumStock <- round(upr2, 2)
  write.csv2(preddata, file = name, row.names = F)
}



# Model checks: all assumptions are well respected
plot(rstandard(lmTest), ylab = "Standardized residuals") # Super-nice looking plot: homoschedasticity respected
resOrdered <- rstandard(lmTest)[order(lamiera$TenDaysOfYear)]
plot(1:nrow(lamiera), resOrdered) # Independence respected
hist(rstandard(lmTest), main="Distribution of standardized residuals", xlab="Standardized residuals")
qqplot(rnorm(length(rstandard(lmTest))) ,rstandard(lmTest), xlab = "Std Normal sample", ylab = "Standardized residuals")
abline(c(0,1), col= "red")# Normality respected
plot(fitted(lmTest), residuals(lmTest, type = "pearson"), xlab = 'Fitted values', ylab = 'Pearson residuals') # mild pattern of pearson residuals, but can be overlooked due to big sample size
# Linearity: respected by construction (i.e. scale of link is always linear)

Business Analytics con retaileR (parte II)

  1. Ecco la seconda parte dell’articolo sul nuovo pacchetto R retaileR. La prima parte potete trovarla qui.

Le Funzioni

closure.opp.cost

closure.opp.cost serve per stimare il costo di opportunità in funzione del tempo di chiusura dei locali di retail. Ossia, una stima del valore monetario perso a seconda dell’orario di chiusura.

I suoi argomenti sono:

Utilizzando il dataset preformattato, ecco un esempio:

closure_opp_cost("16:15:00", X= sales_august)

e l’output:

$Mean
[1] 10.35812

$Median
[1] 10.35

$StDev
[1] 8.392746

$`Late Sales`
[1] 10.30 18.90 6.30 26.70 14.60 10.60 0.05 24.18 10.40 9.60 14.10 14.90 5.10 0.00 0.00 0.00

In questo caso, chiudere il negozio alle 16:15:00 comporterebbe una perdita media di 10.35 €, con una varianza di 8.39 €. Inoltre, dal vettore Late.sales si può vedere chiaramente come 3 giorni non ci siano state vendite dopo le 16:15:00.

 segment_prod_line

Questa funzione è stata creata per l’analisi delle diverse famiglie di prodotti secondo un determinato segmento temporale.

Ossia, se la variabile temporale del dataset input è giornaliera, segment_prod_line produrrà le vendite medie giornaliere per la product line di interesse.

Vediamone il funzionamento:

specifichiamo la linea di prodotti come un vettore di stringhe:

prods &lt;- c("Americano", "Espresso")

e successivamente evochiamo segment_prod_line

segment_prod_line(sales_august, prods)

ecco l’output:

[1] 31.56065

ciò significa che nei giorni che il dataset copre, abbiamo venduto in media 31.5 Americani ed Espressi al giorno.

Mean.items.times

Per ultima, una funzione di analisi temporale delle vendite. mean.items.times calcola le vendite medie (in volume) per ogni ora. Ossia crea un istogramma con l’unità temporale sull’asse delle ascisse e le vendite medie sull’asse delle ordinate.

Ecco il funzionamento (che ricordo, è ottimizzato per un dataset iZettle):

 mean.items.times(sales_august)

e l’output corrispondente:

       Time  Quantity
1  05:00:00  1.000000
2  06:00:00  9.000000
3  07:00:00 23.100000
4  08:00:00 30.695652
5  09:00:00 24.826087
6  10:00:00 17.695652
7  11:00:00 16.000000
8  12:00:00 11.173913
9  13:00:00 10.608696
10 14:00:00  8.590909
11 15:00:00  7.714286
12 16:00:00  6.000000
13 17:00:00  2.000000
14 18:00:00  1.000000

‘Un contributo alla statistica’ di Wislawa Szymborska

No comment: semplicemente una poesia trovata recentemente. A base di statistica ovviamente. Buona lettura.

Un contributo alla statistica di Wislawa Szymborska

Su cento persone:

che ne sanno sempre piu’ degli altri
– cinquantadue;

insicuri a ogni passo
– quasi tutti gli altri;

pronti ad aiutare,
purche’ la cosa non duri molto
– ben quarantanove;

buoni sempre,
perche’ non sanno fare altrimenti
– quattro, be’, forse cinque;

propensi ad ammirare senza invidia
– diciotto;

viventi con la continua paura
di qualcuno o qualcosa
– settantasette;

dotati per la felicita’
– al massimo poco piu’ di venti;

innocui singolarmente,
che imbarbariscono nella folla
– di sicuro piu’ della meta’;

crudeli,
se costretti dalle circostanze
– e’ meglio non saperlo
neppure approssimativamente;

quelli col senno di poi
– non molti di piu’
di quelli col senno di prima;

che dalla vita prendono solo cose
– quaranta,
anche se vorrei sbagliarmi;

ripiegati, dolenti
e senza torcia nel buio
– ottantatre’
prima o poi;

degni di compassione
– novantanove;

mortali
– cento su cento.
Numero al momento invariato.

Business analytics con retaileR (parte I)

Perdonate la lunga assenza: sono stato molto impegnato in questi ultimi mesi. Tuttavia, per farmi perdonare, ho sviluppato una chicca tutta per voi:

un nuovissimo pacchetto R.

L’idea di creare un pacchetto R mi ha sempre affascinato. In fondo, è proprio grazie ai contributi Open Source che questo linguaggio continua ad arricchirsi. Tuttavia, mi è sempre mancato un ambito di applicazione, ossia non ho mai immaginato per cosa costruire il pacchetto. Per fortuna, negli ultimi mesi ho partecipato ad un progetto che mi ha introdotto al mondo del retail (il commercio Business to Consumer per intenderci) ed ai suoi peculiari problemi: da qui è nata l’idea di costruire una serie di funzioni che aiutino l’analisi dei dati in queso settore.

Il pacchetto, sviluppato in R, è chiamato retaileR e potete trovarlo sul mio Github. Nei prossimi mesi cercherò di renderlo disponibile sul CRAN, la repository ufficiale di tutti i pacchetti R.

Che cosa può fare retaileR? Retailer è nato con l’idea di creare un gruppo di funzioni utili per l’analisi dei dati di vendita. In questo modo, il processo decisionale verrebbe reso più semplice e con stime più efficaci. retaileR può aiutare nel decidere l’orario di chiusura di una filiale, organizzare lo staffing individuando gli orari di massimo/minimo volume e molte altre funzioni simili.

In questa serie di articoli illustrerò le funzioni e provvederò degli esempi su come utilizzare questo nuovo pacchetto.


I Dati

Il dataset di base viene considerato quello creato in automatico da iZettle, un sistema POS per il commercio al dettaglio. Qui sotto le prime righe di un dataset illustrativo (che userò negli esempi successivi) e che potrete trovare sempre in GitHub.

Carichiamo il dataset delle vendite fittizie relative ad un mese e visualizziamone le prime righe.

sales_august &amp;lt;- read.csv(“foo/sales_august.csv”)
head(sales_august)

Ecco l’output:

        Date                Time Receipt.number       Name    Variant Unit Quantity Price..GBP. Discount..GBP. Final.price..GBP.
1 2017-08-01 2017-08-20 06:19:00           7359  Americano              NA        1         2.5           0.00              2.50
2 2017-08-01 2017-08-20 06:19:00           7359   Espresso              NA        1         2.0           0.00              2.00
3 2017-08-01 2017-08-20 06:28:00           7360      Mocha              NA        1         3.1           1.55              1.55
4 2017-08-01 2017-08-20 06:43:00           7361 Flat White              NA        2         5.4           0.00              5.40
5 2017-08-01 2017-08-20 06:52:00           7362  Macchiato              NA        1         2.2           0.00              2.20
6 2017-08-01 2017-08-20 06:52:00           7362 Extra Shot Extra shot   NA        1         0.2           0.00              0.20

Le Funzioni

sales.format

questa funzione si propone come l’inizio di ogni analisi con retaileR poiché la sua funzione è formattare correttamente i dati da analizzare. In particolare, essa agisce sulle variabili di tempo (i.e. data e ora) formattandole in maniera uniforme.

Ecco un esempio:

sales_august <- sales.format(sales_sept)
[\code]

I suoi argomenti sono: X, date_var= “Date”, time_var= “Time”, format_time = ‘%H:%M’, dove:

  1. X: indica il data frame contenente i dati da analizzare.
  2. date_var è il nome, come stringa, della variabile contenente la data. Default a “Date”.
  3. time_var idem come date_var ma per il tempo (ossia l’ora). Default a “Time”.
  4. format_time come deve essere formattato il tempo. Per usare retaileR è necessario il formato ‘ora : minuto’ indicato appunto dal default ‘%H:%M’.
as.sales

l’analisi dei dati di vendita crea particolari problemi a livello della variabile temporale. Infatti, essendo una serie temporale, le vendite possono essere raggruppate seguendo tempi diversi: ore, giorni, settimane etc.

Inoltre, non è semplice capire se il di crescita sia positivo o negativo: la stagionalità e i movimenti casuali rendono molto spesso difficile questo compito.

Per evitare questi due problemi, la funzione as.sales crea una nuova classe di oggetti (in S4) che include una semplice regressione lineare e diversi raggruppamenti temporali.

i suoi argomenti sono:

  1. X: indica il data frame contenente i dati da analizzare.
  2. date_var è il nome, come stringa, della variabile contenente la data. Default a “Date”.
  3. time_var idem come date_var ma per il tempo (ossia l’ora). Default a “Time”.
  4. sales_var è il nome, come stringa, della variabile contenete il valore monetario di ogni singola vendita. Defaulta a “Final.price..GBP.”.
  5. format_time come deve essere formattato il tempo. Per usare retaileR è necessario il formato ‘ora : minuto’ indicato appunto dal default ‘%H:%M’.
  6. discount_var, nome colonna contenente gli sconti. Il default è “Discount..GBP.”.

Come output, la funzione produce un oggetto S4 sales, la cui struttura è la seguente:


ecco l'output:

 with 5 slots
  ..@ Daily           :'data.frame':	16 obs. of  2 variables:
  .. ..$ Date             : Date[1:16], format: "2017-08-01" "2017-08-02" "2017-08-03" "2017-08-04" ...
  .. ..$ Final.price..GBP.: num [1:16] 287 268 280 342 337 ...
  ..@ Weekly          :'data.frame':	3 obs. of  2 variables:
  .. ..$ Week             : num [1:3] 31 32 33
  .. ..$ Final.price..GBP.: num [1:3] 1514 2757 1565
  ..@ Monthly         :'data.frame':	1 obs. of  2 variables:
  .. ..$ Month            : num 8
  .. ..$ Final.price..GBP.: num 5836
  ..@ Discount        : num 39.8
  ..@ sales_funct_time:List of 12
  .. ..$ coefficients : Named num [1:2] -201711.1 11.6
  .. .. ..- attr(*, "names")= chr [1:2] "(Intercept)" "Date"
  .. ..$ residuals    : Named num [1:16] 9.79 -21.28 -20.64 29.43 13.04 ...
  .. .. ..- attr(*, "names")= chr [1:16] "1" "2" "3" "4" ...
  .. ..$ effects      : Named num [1:16] -1459 214.3 -18.6 30.4 13.1 ...
  .. .. ..- attr(*, "names")= chr [1:16] "(Intercept)" "Date" "" "" ...
  .. ..$ rank         : int 2
  .. ..$ fitted.values: Named num [1:16] 278 289 301 312 324 ...
  .. .. ..- attr(*, "names")= chr [1:16] "1" "2" "3" "4" ...
  .. ..$ assign       : int [1:2] 0 1
  .. ..$ qr           :List of 5
  .. .. ..$ qr   : num [1:16, 1:2] -4 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:16] "1" "2" "3" "4" ...
  .. .. .. .. ..$ : chr [1:2] "(Intercept)" "Date"
  .. .. .. ..- attr(*, "assign")= int [1:2] 0 1
  .. .. ..$ qraux: num [1:2] 1.25 1.27
  .. .. ..$ pivot: int [1:2] 1 2
  .. .. ..$ tol  : num 1e-07
  .. .. ..$ rank : int 2
  .. .. ..- attr(*, "class")= chr "qr"
  .. ..$ df.residual  : int 14
  .. ..$ xlevels      : Named list()
  .. ..$ call         : language lm(formula = as.formula(paste(sales_var, "~", date_var)), data = day)
  .. ..$ terms        :Classes 'terms', 'formula'  language Final.price..GBP. ~ Date
  .. .. .. ..- attr(*, "variables")= language list(Final.price..GBP., Date)
  .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2] "Final.price..GBP." "Date"
  .. .. .. .. .. ..$ : chr "Date"
  .. .. .. ..- attr(*, "term.labels")= chr "Date"
  .. .. .. ..- attr(*, "order")= int 1
  .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. ..- attr(*, "response")= int 1
  .. .. .. ..- attr(*, ".Environment")=<environment: 0x102165000> 
  .. .. .. ..- attr(*, "predvars")= language list(Final.price..GBP., Date)
  .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "other"
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "Final.price..GBP." "Date"
  .. ..$ model        :'data.frame':	16 obs. of  2 variables:
  .. .. ..$ Final.price..GBP.: num [1:16] 287 268 280 342 337 ...
  .. .. ..$ Date             : Date[1:16], format: "2017-08-01" "2017-08-02" "2017-08-03" "2017-08-04" ...
  .. .. ..- attr(*, "terms")=Classes 'terms', 'formula'  language Final.price..GBP. ~ Date
  .. .. .. .. ..- attr(*, "variables")= language list(Final.price..GBP., Date)
  .. .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. ..$ : chr [1:2] "Final.price..GBP." "Date"
  .. .. .. .. .. .. ..$ : chr "Date"
  .. .. .. .. ..- attr(*, "term.labels")= chr "Date"
  .. .. .. .. ..- attr(*, "order")= int 1
  .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. ..- attr(*, "response")= int 1
  .. .. .. .. ..- attr(*, ".Environment")=<environment: 0x102165000> 
  .. .. .. .. ..- attr(*, "predvars")= language list(Final.price..GBP., Date)
  .. .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "other"
  .. .. .. .. .. ..- attr(*, "names")= chr [1:2] "Final.price..GBP." "Date"
  .. ..- attr(*, "class")= chr "lm"

Come si può osservare, l’oggetto creato da as.sales contiene 5 sub-oggetti:

  •   3 raggruppamenti di tempo: daily, weekly, monthly.
  • Un valore dei discount totali nel dataset.
  • Una oggetto lm di regressione lineare, utile per comprendere il trend.

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!

Quando la Varianza segue la Media: i GLS

Seppur noi umani amiamo etichettare tutti gli eventi e i processi usando categorie ben definite, la realtà è molto più polimorfa. Questa natura delle cose si registra anche in Statistica.

Difatti, abbiamo parlato finora di regressione avendo come modello:

\hat{y} = \beta_{0} + \beta_{1}x + \epsilon ~ N(0,1)

Una delle principali assunzioni di questo modello riguarda gli errori \epsilon che vengono considerati come distribuiti secondo una Normale con media 0 e deviazione standard di 1. In questo caso dunque non esiste alcuna relazione tra la media e la varianza.

Nella vita reale però, lavorando con dati reali, questa assunzione può essere tutto fuorché vera: ci sono molte situazioni in cui tra la media e la varianza esiste una relazione, ossia:

\sigma^{2}_{i} = f(\hat{y_{i}})

In questi casi, ci vengono in aiuto i GLS (Generalized Least Squares). I GLS introducono una forma funzionale esplicita per la varianza degli errori. Le più comuni presentano una relazione positiva fra le due variabili:

  • Relazione di Potenza\hat{y} = \beta_{0} + \beta_{1}x + \epsilon ~ N(0,\sigma^{2}|\hat{y}_{i}|^{2m_{1}})
  • Relazione Esponenziale\hat{y} = \beta_{0} + \beta_{1}x + \epsilon ~ N(0,\sigma^{2}exp(2m_{2}\hat{y_{i}}))

m_{1}m_{2} sono due parametri che devono essere stimati col GLS.

Selezione Automatica di Variabili

In questo articolo vi darò la notizia che ciascuno di voi aspettava. La grandiosità del calcolo matematico unito al potere di calcolo automatico: la selezione automatica di variabili. Ossia, è possibile far selezionare dal computer la combinazione di variabili ottimale per un modello.

L’argomento si innesta su quello dell’articolo precedente, in cui abbiamo illustrato la collinearità.

Esistono tre metodi si selezione automatica di variabili:

  • Eliminazione all’indietro
  • Selezione in avanti
  • Selezione progressiva

Ognuno di questi metodi ha i suoi vantaggi e svantaggi:

Eliminazione all’indietro

Inizia con il creare il modello più complesso possibile, ossia introducendo tutte le variabili disponibili. Ad ogni passo elimina la variabile ‘meno importante’ fino a raggiungere un punto di ottimo usando una certa misura (vedi dopo).

Selezione in avanti

Inizia con il modello più semplice che esiste: solo con l’intercetta. Ad ogni passaggio aggiunge una variabile ‘importante’ fino al raggiungimento di un punto di ottimo.

Selezione progressiva

Combina i due metodi precedenti, cercando ad ogni passo variabili da escludere ed includere. È significativamente più calcolo intensivo.

Se tuttavia vogliamo automatizzare questi algoritmi, ‘l’importanza’ di ogni variabile deve essere misurata in modo oggettivo. Per fare ciò esistono innumerevoli metriche; quella più semplice ed usata si chiama AIC (Aikake Information Criterion). L’AIC è una misura di bontà di predizione penalizzata dal numero di variabili contenute nel modello. Quindi un AIC piccolo segnala un modello migliore. Bisogna tuttavia ricordare che l’AIC è una misura relativa e non assoluta. Formalmente:

AIC = -2 logL + 2p

dove:

  • logL (log-likelihood) è una misura di bontà (ne parleremo in seguito)
  • p è il numero di variabili nel modello statistico.

Un’altra misura molto usata è la BIC (Bayesian Information Criterion), che utilizza un penalità che cambia con la grandezza del campione:

BIC = -2 logL + log(N)p

Dove:

  • è la grandezza totale del campione.

Gli algoritmi di selezione automatica usano quindi queste misura (e altre simili) per dare un rango ad ogni modello (ognuno con una diversa combinazione di variabili) e scegliere quello migliore.

Ovviamente, questi algoritmi automatici sono totalmente indipendenti dal contesto in cui il modello viene creato, i.e. non tengono in considerazione l’ambito scientifico del modello. È qui che il Data Scientist deve usare la sua esperienza ed arte per decidere come agire. L’analista deve dunque capire se prediligere la capacità predittiva (e quindi rischiare di introdurre troppi parametri) o scegliere la capacità descrittiva ed includere solo variabili scientificamente significative.

Largo ad R!

AIC e BIC possono essere calcolati usando queste funzioni, entrambe nel pacchetto {stats}

AIC(object, ...)
BIC(object, ...)

Per la selezione automatica si usa invece questa, sempre nel pacchetto {stats}:

step(object, scope, scale = 0,
direction = c("both", "backward", "forward"),
trace = 1, keep = NULL, steps = 1000, k = 2, ...)

È anche possibile utilizzare la funzione dredge nel pacchetto {MuMIn} che permette la scelta di misure diverse oltre all’AIC.

Ecco un esempio completo di Selezione Progressiva usando il dataset freeny:
Il codice
linearAll = lm(y~., data = freeny)
step.linearALL = step(linearAll, direction = "both")

e l’output

Come vediamo, il modello migliore è con tutte le variabili tranne lag.quarterly.revenue: non includerla produce il più basso AIC possibile.

Collinearità

Una volta costruito un modello lineare col metodo di Maximum Likelihood come si può decidere quali variabili introdurre nel modello? Ovviamente, noi vogliamo introdurre variabili dipendenti che:

  1. abbiano una relazione genuina con la variabile dipendente
  2. aggiungono informazione al modello, considerando le altre indipendenti già presenti nel modello.

Ovviamente, vorremmo allo stesso tempo escludere variabili che forniscono la stessa informazione di quelle già presenti; ossia, vogliamo evitare la collinearità.

Collinearità

Quando variabili collinari sono inserite insieme in un modello, quest’ultimo risulta instabile e otterremo errori standard maggiorati per i parametri. La collinearità può essere individuata tramite i Variance Inflation Factors (VIFs).

I VIFs si calcolano costruendo un modello lineare tra ogni covariata e tutte le altre covariate. Ad esempio, per la prima variable indipendente

x_{1i} = \delta_{0} + \delta_{1}x_{2i} + \cdots + \delta_{p}x_{pi} + u_{i}

dove u_{i} è l’errore.

Un modello di questo tipo è costruito per ogni indipendente e prendendo l’R2 si calcola il VIF:

VIF_{j} = \frac{1}{1-R^{2}_{j}}

per le covariate j = 1, . . ., p.

Grandi VIFs indicano collinearità. Tuttavia non ci sono regole ben precise sul ‘grandi’: la regola più comune è di considerarli problematici quando VIF > 4, perché l’intervallo di confidenza (IC) del parametro j sarà grande il doppio del normale (ossia \sqrt{VIF} indica quanto sarà dilatato l’IC)

Facciamo un esempio in R:

LinearAll = lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data = iris) # Creiamo il modello lineare
library(car) # Carichiamo la libreria per la funzione 'vif'
vif(LinearAll)

Ecco l’output:

Sepal.Width Petal.Length Petal.Width
1.270815    15.097572   14.234335

Come possiamo vedere le variabili ‘Petal.Length’ e ‘Petal.Width’ sono altamente collineari: se le togliessimo otterremo una maggiore stabilità del modello

Un’altra prova del fatto che che il modello sia affetto da collinearità sono parametri il cui valore cambia sostanzialmente quando nuove variabili sono aggiunte al modello.

Come risolvere la collinearità? La soluzione più semplice è di rimuovere le variabili collinari: molto spesso questo accade automaticamente quando si usano algoritmi di selezione automatica.

Post-Verità e la casualità dell’esistenza

O anche: ‘La Statistica è morta. Lunga vita alla Statistica!’


I recenti sviluppi politici mondiali hanno sicuramente confuso i più riguardo cosa sia cambiato nella mente dell’elettorato, sul perché idee apparentemente irrazionali1 e populiste abbiano avuto il sopravvento. L’evento più lampante (se non folkloristico) è stata l’esclamazione di Gove, un politico sostenitore della Brexit:

Britain has had enough of experts

Questa frase, oltre ad avere un significato filosofico profondo, ha anche dei risvolti sinificativi a livello statistico, perché riflette la sfiducia del cittadino contemporaneo nella misurazione quantitativa dello realtà, ossia la Statistica.

Un recente articolo apparso sul Guardian ha per l’appunto affrontato questo tema, che risulta essere epocale, tanto da aver generato il neologismo di post-verità.

Uno studio di Marketplace ha rilevato come il 47.5% dell’elettorato Trump e il 25% del totale “non crede ai dati economici ufficiali pubblicati dal governo federale”.

Risultati studio Marketplace sull’elettorato americano

Allo stesso modo, in un altro studio YouGov/Uni Cambridge, il 55% dei rispondenti crede che “il governo stia nascondendo il numero reale di immigrati viventi nel paese”.

Pare dunque ovvio che il pubblico occidentale non creda più all’oggettività della misurazione e al metodo scientifico della Statistica. Non solo, sembra anche, come risulta dalle parole di Gove, che tale oggettività risulti arrogante e debba essere combattuta; ecco dunque che gli unici argomenti che contano sono solo soggettivi e altamente personalizzati.

Un esempio? Un report del think-tank BritishFuture ha dimostrato come le persone siano emotivamente colpite da storie personali e struggenti di immigrati; al contrario, i “freddi” dati generano l’effetto diametralmente opposto, soprattutto se dimostrano l’influenza positiva degli immigrati sull’economia. La motivazione? Le persone assumono automaticamente che siano dati contraffatti.

Si potrebbe discutere per ore sulle motivazioni profonde di questa ostilità verso i dati statistici ufficiali (una mia spiegazione personale e ‘a pelle’ è che il bias di conferma sia più forte di quanto pensiamo), ma forse questa non è la sede più adatta.

Tuttavia, voglio sottoporci ad un esperimento mentale riguardo le implicazioni statistiche di questo fenomeno: e se avessero ragione? Ossia, se i dati fossero davvero non rappresentativi della realtà?

Questa domanda fa riflettere in maniera profonda sulla natura dei dati statistici: essi sono infatti nati per riassumere la stocasticità del mondo in maniera oggettiva ed inequivocabile. Esistono proprio per demolire la diversità e riassumerla in una piccola panoplia di numeri in modo tale da ‘rendere’ la realtà più semplice ed interpretabile dalle nostre menti limitate. Il problema però, è che seppur la nostra rappresentazione del mondo è ora semplice, il mondo rimane complesso, qualsiasi modello tu voglia usare. Esempio lampante in statistica economica: il GDP misura davvero la produzione di un paese? A che geografia è meglio produrlo? Nazionale? Regionale? Che cosa vogliamo contarci dentro: produzione metallurgica? Prostituzione? E i lavori domestici?

Ecco dunque l’insegnamento centrale della Statistica, che a mio avviso dovrebbe essere insegnato prima ancora dei numeri aleatori:

la Statistica è un compromesso tra il comprendere e il descrivere.

Più la descrizione è dettagliata, meno è comprensibile, e vice versa. Non a caso esistono misure come R2 o la Devianza, proprio per aiutarci a scegliere un livello di questo compromesso. Ma ancora: sono di nuovo misure, usate per misurare altre misure. Sta al Data Scientist, allo Statistico l’onere di scegliere, di valutare tutte queste misure, perché da soli i numeri non vogliono dire nulla: la magia sta tutta qui, nell’esperienza dell’analista che usa una mente umana per cercar di comprendere la stocasticità dell’esistenza.

Resta solo da vedere che strada prenderanno le nuove tendenze Statistiche come la Machine Learning e la Data Science, che creano delle vere e proprie macchine automatiche: la seconda vita della Statistica; dai semplice numeri, alle macchine.

E intanto, il mondo resta comunque casuale, in barba a tutti noi.

Note

1. Seppur io faccia fatica a non esprimermi a livello politico, cerco di resistere alla brama.