Esercitazione di Social Media Management A.A. 2016­2017Regressione Lineare

Antonino Furnari http://www.dmi.unict.it/~furnari/ ­ furnari@dmi.unict.it Prof. Giovanni Maria Farinella http://www.dmi.unict.it/~gfarinella/ ­ gfarinella@dmi.unict.it

Si consiglia la consultazione online

In questa esercitazione impareremo come:

Analizzare un insieme di osservazione mediante Principal Component Analysis (PCA);Costruire un modello di regressione lineare per predire la “memorability” di una immagine apartire da features di tipo Bag Of Visual Words;Stimare l’accuratezza del modello mediante errore assoluto medio (Mean Absolute Error ­ MAE) eindice di correlazione per ranghi di Spearman.

1. ConvenzioniLo scopo di questa esercitazione è guidare il lettore attraverso la comprensione di concetti chiave e la loroimplementazione pratica. Per facilitare la comprensione degli argomenti proposti, durante l'esercitazione,viene richiesto al lettore di rispondere ad alcune domande e risolvere alcuni esercizi pratici. Tali richiestesono indicate dai seguenti simboli:

Questo simbolo è presente laddove viene richiestoal lettore di rispondere a una domanda.

Questo simbolo è presente laddove viene richiesto alettore di svolgere un esercizio.

Questo simbolo indica una nota al lettore.

Questa icona contrassegna i contenuti nonessenziali ai fini dell'esercitazione. L'esercitazionepuò essere svolta in modalità "fast track" saltando icontenuti contrassegnati da questo simbolo. L'iconaè generalmente riportata in alto a sinistra nelparagrafo interessato.

Questa icona indica un collegamento a una pagina di approfondimento. Si trattadi materiale non indispensabile ai fini dell'esercitazione ma che può comunquetornare utile. L'icona è generalmente riportata in alto a destra nel paragrafointeressato, cliccandovi sopra è possibile accedere al contenuto diapprofondimento. Gli approfondimenti segnalati da questa icona sono ancheriportati tra le referenze.

Durante l'esercitazione, verrà richiesto di sviluppare alcune funzioni necessarie a svolgere i passisuccessivi. Le soluzioni non sono riportate in questa esercitazione

Le funzioni richieste si trovano nel modulo solution che verrà fornito al completamento dell'esercitazione.

In [1]:

#importiamo tutte le funzioni del package solution che ci serviranno in seguito #Nota: questa istruzione è qui solo per "motivi tecnici" #Ai fini di questa esercitazione, non è richiesta l'esecuzione di questa istruzione #che non andrebbe a buon fine vista l'assenza del modulo "solution" from solution import *

2. RequisitiAmbiente di lavoro basato su python 2.7 (consigliato IPython + IDE). L'ambiente di lavoro puòessere impostato installando una delle principali distribuzioni Python:

Anaconda (https://www.continuum.io/downloads);Canopy (https://www.enthought.com/products/canopy/).

Librerie richieste:scipy (https://www.scipy.org/);numpy (http://www.numpy.org/);matplotlib (http://matplotlib.org/);scikit­image (http://scikit­image.org/).

È inoltre richiesto aver seguito le seguenti esercitazioni:

Gestione di dataset di immagini: http://www.dmi.unict.it/~furnari/teaching/SMM1617/lab0/;Bag of Visual Word Models: http://www.dmi.unict.it/~furnari/teaching/SMM1617/lab1/;

In particolare, il file dataset.py (scaricabile da qui) deve essere posizionato nella directory di lavoro.Utilizzeremo anche le funzioni sviluppate nell'esercitazione su Bag of Visual Word Models, che sonoriportate nel file bovw.py (scaricabile da qui). Posizionare anche questo file nella directory di lavoro.

3. DatasetIn questa esercitazione utilizzeremo un dataset di immagini acquisite nel contesto di uno studio su“Image Memorability” (vedi referenza [1]). Il dataset contiene 2222 immagini di risoluzione 256×256pixels. Le immagini sono corredate da una serie di valori stimati mediante un gioco di memoria visuale alquale hanno partecipato diversi soggetti utilizzando la piattaforma Amazon Mechanical Turk.

Tra i valori forniti, noi considereremo gli “hits”, ovvero il numero di volte in cui ogni immagine è stataeffettivamente riconosciuta come “già vista” quando ripresentata ai partecipanti. Il dataset è disponibileall’indirizzo

http://web.mit.edu/phillipi/Public/WhatMakesAnImageMemorable/

Tuttavia, considerate le dimensioni dell’archivio (1.7 GB), per questa esercitazione utilizzeremo unaversione “più leggera” del dataset in cui tutte le immagini sono state compresse in formato JPEG. Il datasetè già stato diviso in training e test set dopo un riordinamento casuale delle immagini. Per risparmiaretempo, ai fini di questa esercitazione, vengono inoltre fornite le rappresentazioni estratte dalle immagine apartire da un modello di tipo Bag of Visual Words.

Scarichiamo il dataset da http://yoda.dmi.unict.it/memorability_python.zip e scompattiamo l’archivio nellanostra directory di lavoro. All’interno della cartella "memorability" troveremo:

Una cartella "training_set" contenente una cartella "Images" con 1555 immagini in formato JPEG(immagini di training). Il nome di ogni immagine rappresenta un indice numerico che definisce unordinamento univoco;Una cartella "test_set" contenente una cartella "Images" con 667 immagini in formato JPEG(immagini di test). Il nome di ogni immagine rappresenta un indice numerico che definisce unordinamento univoco;Il file “training_hits.npy”, che contiene il numero “hits” per ogni immagine di training. I valori sonoordinati secondo la numerazione definita dai nomi dei file JPG;Il file “test_hits.npy”, che contiene il numero “hits” per ogni immagine di test. I valori sono ordinatisecondo la numerazione definita dai nomi dei file JPG;Il file “features.pkl”, che contiene le features di tipo Bag Of Visual Words pre­computate per ogniimmagine come specificato di seguito.

Iniziamo costruendo due oggetti di tipo Dataset a partire dalle immagini di training e test appena estrattedall'archivio:

In [2]:

from dataset import Dataset training_set = Dataset('memorability/training_set') test_set = Dataset('memorability/test_set')

La classe Dataset è stata progettata per gestire dataset di immagini appartenentia più classi. In questa esercitazione, affronteremo un problema di regressione enon di classificazione. Pertanto, il dataset conterrà una unica "pseudo­classe"chiamata "Images".

Verifichiamo che il dataset sia stato correttamente caricato e visualizziamo alcune immagini:

In [3]:

print "Number of images in the training set:0".format(training_set.getLength()) print "Number of images in the test set:0".format(test_set.getLength()) print "Total number of images:0".format(training_set.getLength()+test_set.getLength()) print "List of classes ", training_set.getClasses() training_set.showImage('Images',18) print "Shown image path: %s" % training_set.getImagePath('Images',18)

Carichiamo dunque il numero di "hits" per ogni immagine del training set:

In [4]:

import numpy as np training_hits=np.load('memorability/training_hits.npy') test_hits=np.load('memorability/test_hits.npy') print training_hits.shapeprint test_hits.shape

Possiamo dunque ispezionare il numero di hits relativo a una data immagine:

Number of images in the training set:1555 Number of images in the test set:667 Total number of images:2222 List of classes ['Images']

Shown image path: memorability/training_set\Images\0019.jpg

(1555L,) (667L,)

In [5]:

from skimage import io as sio from matplotlib import pyplot as plt im=sio.imread(training_set.getImagePath('Images',23)) plt.figure() plt.imshow(im) plt.title("Hits: 0".format(training_hits[23])) plt.show()

Esercizio 3.1 Scrivere il codice per visualizzare tutte le immagini del training set con i rispettivinumeri di hits in sequenza. Suggerimento: utilizzare la funzione `plt.pause(n)` tradue esecuzioni consecutive di `imshow` in modo da visualizzare l'immagine per`n` secondi.

Domanda 3.1 Notate una correlazione tra contenuto visivo dell’immagine e numero di hits?

Esercizio 3.2 Visualizzare le immagini con il numero più basso e più alto di hits. Suggerimento:vedere la documentazione della funzione di numpy argmax.

4. RappresentazioneAdesso eseguiamo tutti i passi necessari per creare un modello di tipo Bag Of Visual Words e estrarre lefeature da tutto il dataset. Ci avvarremo delle funzioni sviluppate nell'esercitazione precedente.

La costruzione di un modello di tipo Bag Of Visual Words e la rappresentazione ditutte le immagini del training e test set possono essere dispendiose in termini ditempo. Pertanto, ai fini di questa esercitazione, sono state fornite le feature pre­computate per ogni immagine di training e test. E' dunque possibile saltare ilpunto 4.1 (o leggerlo velocemente) per passare direttamente al punto 4.2. Siconsiglia di ritornare sul punto 4.1 a casa, una volta ultimata l'esercitazione.

4.1 Rappresentazione delle immagini e salvataggio su file (dafare a casa)

Per creare il modello di tipo Bag Of Visual Words, utilizzeremo alcune delle funzioni sviluppatenell'esercitazione su Bag Of Visual Words model. Per ottenere una rappresentazione più accurata,utilizzeremo step=4 (valore di default) nel processo di estrazione delle feature locali DAISY. Le funzioni,modificate in modo da utlizzare il corretto valore di step sono definite all'interno del file bovw.py(scaricabile da qui), che va posizionato nella directory di lavoro.

from bovw import extract_features training_local_features = extract_features(training_set) from bovw import describe_dataset from sklearn.cluster import MiniBatchKMeans as KMeans kmeans = KMeans(500) kmeans.fit(training_local_features) training_features,_,_=describe_dataset(training_set,kmeans) test_features,_,_=describe_dataset(test_set,kmeans)

Normalizziamo le feature mediante norma L2:

from sklearn.preprocessing import Normalizer norm = Normalizer(norm='l2') training_features=norm.transform(training_features) test_features=norm.transform(test_features)

Quando si affrontano problemi di regressione, è comune chiamare le feature e i target . Pertantoponiamo:

X_training = training_features y_training = training_targets["Images"] X_test = test_features y_test = test_targets["Images"]

Salviamo dunque le feature su file:

import cPickle with open("memorability/features.pkl","wb") as f: cPickle.dump( "X_training" : X_training, "y_training" : y_training, "X_test" : X_test, "y_test" : y_test ,f)

4.2 Caricamento delle feature da file (se si è saltato il passo4.1)Dato che l’operazione può essere molto dispendiosa in termini di tempo computazionale, carichiamo ilmodello e le feature pre­computate per questa esercitazione:

X y

In [6]:

import cPickle with open("memorability/features.pkl") as f: data = cPickle.load(f) X_training = data["X_training"] y_training = data["y_training"] X_test = data["X_test"] y_test = data["y_test"]

Domanda 4.1 Che dimensioni hanno le matrici X_training e X_test? Questo tipo di matrice ha unnome particolare in machine learning, quale?

5. Principal Component AnalysisQuando si utilizza l’equazione normale per trovare i parametri ottimali di un regressore lineare, sipossono avere problemi relativi alla non­invertibilità di , dove è la matrice delleosservazioni. Ciò è generalmente dovuto alla presenza di “ridondanza” tra le feature (diverse parole visualiche codificano pattern visivi molto simili, nel nostro caso) o all’eccessivo numero di feature rispetto alnumero di osservazioni. Per prevenire questo tipo di problema, trasformeremo i nostri dati utilizzando unatecnica chiamata Principal Component Analysis (PCA).

Data una matrice di osservazioni (con numero di osservazioni e numero di feature), laPCA trova una matrice tale che:

dove è la matrice delle osservazioni centrata nello zero ( , con valore medio di neltraining set), è la pseudo­inversa di e è una matrice contenente osservazioni didimensionalità generalmente ridotta .

In questa sede, non scenderemo nei dettagli matematici della PCA (per i quali si è rimandati alla relativapagina di wikipedia. Ci basta solo sapere che la nuova matrice delle osservazioni , ha le seguenticaratteristiche:

Le features sono ordinate dalla “più importante” alla “meno importante”. Ciò vuol dire che lamaggior parte dell’informazione è contenuta nelle prime features;Le features sono generalmente “decorrelate” e quindi sono meno inclini a presentare dipendenzalineare.

Le medie e la matrice di proiezione dipendono dai dati di training che abbiamo utilizzato per stimarli.Si suppone che, se la trasformazione PCA definita da questi valori è "buona" per i dati di training, essasarà altrettanto "buona" per quelli di test. Pertanto, una volta stimate le medie e la matrice , èpossibile applicare la stessa trasformazione PCA al training e al test set in modo da avere dati dallecaratteristiche analoghe.

Il calcolo dei valori di e e l'applicazione di una data trasformazione PCA su nuovi dati, può essereottenuta mediante l'oggetto PCA di sklearn. In particolare, possiamo "allenare" il modello PCA mediante leseguenti righe di codice:

In [7]:

from sklearn.decomposition import PCA as PCA pca = PCA() pca.fit(X_training);

Il metodo fit ha stimato le medie e la matrice a partire dai dati di training passati in input. E' possibileaccedere ai valori stimati come segue:

( XXT )−1 X

X[n × m] n m

C[m × p]

W = (X¯ ¯¯

C T )−1

X¯ ¯¯ = −X

¯ ¯¯i Xi μi μi Xi

(C T )−1 C T W [n × p] n

p ≤ m

W

μi C

μi C

μi C

μi C

In [8]:

mu = pca.mean_ C = pca.components_ print mu.shape print C.shape

A questo punto, per applicare la trasformazione PCA imparata, potremmo utilizzare la formula vista prima:

che si traduce in python come segue:

In [9]:

W=np.dot((X_training‐mu),np.linalg.inv(C)) print W.shape

Tuttavia, mediante la libreria scikit­learn è sufficiente scrivere:

In [10]:

X_training_pca = pca.transform(X_training) X_test_pca = pca.transform(X_test) print X_training.shape, X_training_pca.shape print X_test.shape, X_test_pca.shape

In questo caso le dimensioni di X_training_pca e X_test_pca sono rimaste invariate (500 feature). Tuttavia,le feature sono adesso state ordinate in maniera "più significativa".

Per capire in che modo sono stati trasformati i dati, proviamo a plottare in uno spazio 2D i campionipresenti in X_training e X_training_pca selezionando solo le prime due componenti:

W = (X¯ ¯¯

C T )−1

(500L,) (500L, 500L)

(1555L, 500L)

(1555L, 500L) (1555L, 500L) (667L, 500L) (667L, 500L)

In [11]:

from matplotlib import pyplot as plt plt.figure(figsize=(10,4)) plt.subplot(121) plt.plot(X_training[:,0],X_training[:,1],'bx') plt.subplot(122) plt.plot(X_training_pca[:,0],X_training_pca[:,1],'rx') plt.show()

Domanda 5.1 Che differenze notate tra i due plot? Quale insieme appare “centrato” rispetto allozero? Quale insieme sembra contenere meno ridondanza?

In pratica, questa nuova rappresentazione delle osservazioni ci permette di:

Avere meno problemi quando si utilizza la regressione in forma normale;Avere un’idea del modo in cui i nostri dati si dispongono nello spazio selezionando solo le primedue componenti.

6. Costruire un regressore lineare multivariato perpredire la “memorability” di una immagine

Alleneremo un modello di regressione lineare che, presa in input la rappresentazione di unaimmagine ne stimi il numero di hits (ovvero la "memorability") mediante ilmodello:

Nell'equazione riportata sopra, è la variabile dipendente (numero di hits), sono i predittori (le 500 feature estratte mediante BOVW e trasformate secondo PCA) è l’intercetta, e

sono i regressori. Per poter utilizzare il modello di regressione riportato sopra, dobbiamoprima "allenarlo", cioè stimare i parametri (regressori) a partire dalle immagini del training set.

Per allenare il modello, utilizzeremo l'oggetto LinearRegression della librearia scikit­learn che ha unainterfaccia analoga a quella di PCA:

In [12]:

from sklearn.linear_model import LinearRegression lr = LinearRegression() lr.fit(X_training_pca,y_training);

E' possibile accedere ai parametri stimati come segue:

In [13]:

theta_0 = lr.intercept_ theta_1_500 = lr.coef_ print theta_0 print theta_1_500.shape

Domanda 6.1 Perché il parametro viene chiamato "intercept_" all'interno del modello diregressione lineare?

Se vogliamo rappresentare graficamente come il numero di hits varia al variare delle osservazioni e comesi comporta il nostro regressore lineare, selezioneremo un solo predittore, ignorando gli altri. Dal momentoche la PCA ha “ordinato” le feature per importanza, considereremo la prima feature. Il modello si riduce a:

Visualizziamo innanzitutto come la variabile dipendente varia al variare del predittore :

X = ( , , … , )X1 X2 Xn y

y = + + + … +θ0 X1θ1 X2θ2 X500θ500

y , … ,X1 X500 500θ0

, … ,θ1 θ500

, … ,θ0 θ500

θ0

y = +θ0 X1θ1

y X1

52.3009646302 (500L,)

In [14]:

plt.figure() X1 = X_training_pca[:,0] plt.plot(X1,y_training,'x') plt.show()

A questo punto possiamo visualizzare la retta di regressione definita dal nostro modello linearesemplificato:

In [15]:

plt.figure() plt.plot(X1,y_training,'x') theta_1=theta_1_500[0] x=np.arange(‐0.6,0.8,0.1) plt.plot(x,theta_0+x*theta_1,'r') plt.show()

Dato che stiamo effettuando una regressione multivariata, questo plot “cattura” cosa succederebbe se ivalori delle feature superiori alla prima fossero nulli.

Per avere un’idea di come la regressione “si comporta” all’aumentare del numero di feature, passiamodallo spazio bidimensionale a quello tridimensionale considerando i primi due predittori e . Il nostromodello lineare diventa:

Visualizziamo i dati in tre dimensioni, mediante i comandi:

X1 X2

Y = + +θ0 X1θ1 X2θ2

In [16]:

from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() plt.subplot(111, projection='3d') X1 = X_training_pca[:,0] X2 = X_training_pca[:,1] plt.plot(X1,X2,y_training,'o') plt.show()

Il plot appena ottenuto è un plot tridimensionale. Pertanto è possibile "esplorarlo"cambiandone il punto di vista. Per farlo, provate a cliccare su un punto qualsiasidel grafico e trascinare con il mouse. Si noti che, nel caso in cui si utilizza jupyter,sarà generata una immagine statica. Per ottenere una versione esplorabile delgrafico si digiti: ```python %matplotlib qt ``` prima di lanciare il plot.

Per poter effettuare i plot in 3 dimensioni, abbiamo inserito un subplot all'interno della figura appena creatae utilizzato la keyword projection='3d' per specificare che il plot sarebbe stato in 3 dimensioni. Permaggiori informazioni sui plot 3D è possibile consultare la documentazione di matplotlib.

Il modello di regressione lineare in 3 dimensioni non è più una retta ma un piano individuatodall'equazione:

Dove e sono due variabili indipendenti, mentre è la variabile dipendente. Possiamo visualizzare ilpiano come segue:

z = + x + yθ0 θ1 θ2

x y z

In [17]:

fig = plt.figure() plt.subplot(111, projection='3d') theta_2 = theta_1_500[1] x_range= np.arange(‐1,1,0.1) y_range= np.arange(‐1,1,0.1) x,y=np.meshgrid(x_range,y_range) plt.plot(X1,X2,y_training,'x') plt.gca().plot_surface(x,y,theta_0+theta_1*x+theta_2*y,shade=False,color='r') plt.show()

La funzione “meshgrid” permette di creare una griglia bidimensionale di valori a partire da due vettoriunidimensionali.

In [18]:

x=[1,2,3,4,5] y=[1,2,3,4,5,6,7,8] [x_grid,y_grid]=np.meshgrid(x,y) print "x:",x print "y:",y print "x_grid:\n",x_grid print "y_grid:\n",y_grid

La funzione “plot_surface” permette di plottare una superficie nel mondo 3D (un piano nel nostro caso).

Domanda 6.2 la retta e il piano di regressione sono “ragionevoli” rispetto ai dati considerati?Quanto errore pensate che il regressore stia compiendo?

Esercizio 6.1 Visualizzare i due grafici come quelli di sopra anche per il test set. Che risultatiabbiamo? La retta e il piano di regressione “si comportano” in maniera similerispetto ai nuovi dati?

x: [1, 2, 3, 4, 5] y: [1, 2, 3, 4, 5, 6, 7, 8] x_grid: [[1 2 3 4 5] [1 2 3 4 5] [1 2 3 4 5] [1 2 3 4 5] [1 2 3 4 5] [1 2 3 4 5] [1 2 3 4 5] [1 2 3 4 5]] y_grid: [[1 1 1 1 1] [2 2 2 2 2] [3 3 3 3 3] [4 4 4 4 4] [5 5 5 5 5] [6 6 6 6 6] [7 7 7 7 7] [8 8 8 8 8]]

In [19]:

fig = plt.figure() plt.subplot(111, projection='3d') X1=X_test_pca[:,0] X2=X_test_pca[:,1] theta_2 = theta_1_500[1] x_range= np.arange(‐1,1,0.1) y_range= np.arange(‐1,1,0.1) x,y=np.meshgrid(x_range,y_range) plt.plot(X1,X2,y_test,'x') plt.gca().plot_surface(x,y,theta_0+theta_1*x+theta_2*y,shade=False,color='r') plt.show()

A questo punto possiamo ottenere le predizioni relative a training e test set mediante il metodo “predict”:

In [20]:

training_predictions = lr.predict(X_training_pca) test_predictions = lr.predict(X_test_pca)

Esercizio 6.2 Scrivere una funzione predict_and_show che prenda in input:

un oggetto di tipo Dataset;le rappresentazioni delle immagini del dataset (già trasformatemediante PCA);il modello di regressione imparato;i valori target relativi alle immgini del dataset (ad esempio `y_training` o`y_test`);un indice numerico;

e mostri l'immagine del dataset corrispondente all'indice specificato, con ilnumero di hits predetto e il relativo valore di ground truth. La funzione dovrà dareun risultato simile al seguente:

X

In [21]:

predict_and_show(training_set,X_training_pca,lr,y_test,1) predict_and_show(test_set,X_test_pca,lr,y_test,18)

Domanda 6.3 Il regressore lineare che abbiamo allenato funziona? Come potremmo valutarnela bontà?

7.2 Mean Absolute Error (MAE)Il Mean Absolute Error (MAE) è l’errore assoluto medio che compiamo nelle predizioni. In pratica, date levariabili target e le predizioni restituite dal nostro modello, il MAE è definito come:

Esercizio 7.1 Calcolare il MAE relativo a training e test set. Quale valore è più basso?

Domanda 7.1 Considerando che il MAE rappresenta l’errore medio in valore assoluto checompiamo quando cerchiamo di predire il numero di “hits” che una immaginetotalizzerebbe nel memory game, l’errore che commettiamo è alto o basso?

7.3 Indice di correlazione per ranghi di SpearmanL’indice di correlazione per ranghi di Spearman indica quanto due variabili siano dipendenti. Si trattadi un numero compreso tra ­1 (massima anticorrelazione) e 1 (massima correlazione). Secondoquesto indice, due variabili sono correlate positivamente (rips. negativamente) quando una può essereespressa come una funzione monotona crescente (risp. decrescente) dell’altra. Più la monotonia è“perfetta” più il numero si avvicina a 1 o ­1, come mostrato negli esempi di seguito:

L'indice di correlazione di spearman può essere calcolato mediante la funzione spearmanr di scipy (siveda la documentazione):

yi y i

MAE( , ) = | − |yi y i

1

n∑i=1

n

yi y i