# Classificazione di vini attraverso gaussiane unidimensionali¶

Il dataset Wine sarà il nostro esempio ricorrente per la discussione dell' approccio generativo alla classificazione.

Il dataset (costruito nel 1991 dall'Università di Genova) può essere scaricato dal repository UCI (https://archive.ics.uci.edu/ml/datasets/wine). Contiene 178 osservazioni, ognuna corrispondente ad una bottiglia di vino:

- Le feature (`x`): un vettore 13-dimensionale consistente di caratteristiche chimiche e visuali della bottiglia di vino
- Le etichette (`y`): la cantina di provenienza della bottiglia (1,2,3)

Prima di continuare, assicurarsi che il dataset (`wine.data.txt`) sia presente in una sottocartella di nome `data`.

## 1. Caricamento del dataset¶

Iniziamo caricando i pacchetti Python che useremo.

In [1]:

    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    # Un modulo utile per la gestione di distribuzioni gaussiane
    from scipy.stats import norm, multivariate_normal
    # Moduli per la visualizzazione di grafici interattivi
    import ipywidgets as widgets
    from IPython.display import display
    from ipywidgets import interact, interactive, fixed, interact_manual, IntSlider

Ora carichiamo il dataset Wine. Ci sono 178 osservazioni, ciascuna con 13 feature ed un'etichetta (1,2,3). Le divideremo in un training set di 130 punti ed un test set di 48 punti.

In [2]:

    data = np.loadtxt('data/wine.data.txt', delimiter=',')
    # Nomi delle feature
    featurenames = ['Alcool', 'Acido malico', 'Ceneri', 'Alcalinità delle ceneri', 'Magnesio', 'Fenoli totali',
                   'Flavonoidi', 'Fenoli non flavonoidi', 'Proantocianina', 'Intensità di colore', 'Tonalità',
                  'OD280/OD315', 'Prolina']

Fissiamo una particulare permutazione pseudocasuale dei dati, ed usiamola per effettuare la partizione in training e test. Vogliamo creare 4 array:

- `trainx`: 130x13, le osservazioni di training
- `trainy`: 130x1, le etichette delle osservazioni di training
- `testx`: 48x13, le osservazioni di test
- `testy`: 48x1, le etichette delle osservazioni di test

In [3]:

    # Suddividiamo le 178 istanze in un training set (trainx, trainy) di taglia 130 e un test set (testx, testy) di taglia 48
    # Separiamo inoltre i dati di input (colonne 1-13) dalle etichette (colonna 0)
    np.random.seed(0) # inizializzazione del generatore pseudo-casuale (per la riproducibilità degli esperimenti)
    perm = np.random.permutation(178)
    trainx = data[perm[0:130],1:14]
    trainy = data[perm[0:130],0]
    testx = data[perm[130:178], 1:14]
    testy = data[perm[130:178],0]

Verifichiamo quanti dati di training ci sono in ogni classe.

In [4]:

    sum(trainy==1), sum(trainy==2), sum(trainy==3)

Out[4]:

    (43, 54, 33)

### Esercizio rapido 1¶

Quante istanze di test sono presenti in ciascuna classe?

In [5]:

    # modificare questa cella

Cliccare qui per la soluzione.

## 2. Costruzione della distribuzione di una singola feature da ognuna delle cantine¶

Focalizziamoci su una singola feature: 'Alcool'. Questa è la prima feature, quella di indice 0. Il seguente grafico è un'istogramma dei valori di questa feature nella classe 1 (cantina 1), insieme al fit gaussiano della distribuzione.

Come si può generare un grafico di questo tipo?

La seguente funzione, density_plot, lo fa per ogni feature ed etichetta. La prima riga aggiunge una componente interattiva che ci permette di scegliere questi parametri utilizzando dei controlli.

Provatelo, e poi, osservate il codice attentamente per capire esattamente ciò che fa ciascuna riga della funzione.

In [6]:

    @interact_manual( feature=IntSlider(0,0,12), label=IntSlider(1,1,3))
    def density_plot(feature, label):
        plt.hist(trainx[trainy==label,feature], density=True)
        #
        mu = np.mean(trainx[trainy==label,feature]) # media
        var = np.var(trainx[trainy==label,feature]) # varianza
        std = np.sqrt(var) # deviazione standard
        #
        x_axis = np.linspace(mu - 3*std, mu + 3*std, 1000)
        plt.plot(x_axis, norm.pdf(x_axis,mu,std), 'r', lw=2)
        plt.title("Cantina "+str(label) )
        plt.xlabel(featurenames[feature], fontsize=14, color='red')
        plt.ylabel('Densità', fontsize=14, color='red')
        plt.show()

### Esercizio rapido 2¶

Domanda: per quale feature (indice da 0 a 12) la distribuzione dei valori (nel training set) della Cantina 1 ha la deviazione standard minima?

In [7]:

    # modificate questa cella

Cliccare qui per la soluzione.

## 3. Fit di una gaussiana per ogni classe¶

Definiamo una funzione che effettua il fit di un modello generativo gaussiano alle tre classi, restringendosi ad un'unica feature.

In [8]:

    # La funzione assume che l'etichetta y vari nell'insieme {1,2,3}
    def fit_generative_model(x,y,feature):
        k = 3 # numero di classi
        mu = np.zeros(k+1) # lista delle medie
        var = np.zeros(k+1) # lista delle varianze
        pi = np.zeros(k+1) # vettore dei pesi delle classi
        for label in range(1,k+1):
            indices = (y==label)
            mu[label] = np.mean(x[indices,feature])
            var[label] = np.var(x[indices,feature])
            pi[label] = float(sum(indices))/float(len(y))
        return mu, var, pi

Chiamiamo questa funzione sulla feature di indice 0 (Alcool). Quali sono i pesi delle classi?

In [9]:

    feature = 0 # 'alcool'
    mu, var, pi = fit_generative_model(trainx, trainy, feature)
    print(pi[1:])
    print(mu[1:])
    print(var[1:])

    [0.33076923 0.41538462 0.25384615]
    [13.78534884 12.31092593 13.15969697]
    [0.23325279 0.2819047  0.2851787 ]

Grafichiamo la distribuzione Gaussiana risultante per ciascuna delle tre classi.

In [10]:

    @interact_manual( feature=IntSlider(0,0,12) )
    def show_densities(feature):
        mu, var, pi = fit_generative_model(trainx, trainy, feature)
        colors = ['r', 'k', 'g']
        for label in range(1,4):
            m = mu[label]
            s = np.sqrt(var[label])
            x_axis = np.linspace(m - 3*s, m+3*s, 1000)
            plt.plot(x_axis, norm.pdf(x_axis,m,s), colors[label-1], label="Cantina " + str(label))
        plt.xlabel(featurenames[feature], fontsize=14, color='red')
        plt.ylabel('Densità', fontsize=14, color='red')
        plt.legend()
        plt.show()

### Esercizio rapido 3¶

Usando il grafico interattivo soprastante, rispondere alle seguenti domande:

- Per quale feature (0-12) le distribuzioni della classe 1 e 3 si sovrappongono di più?
- Per quale feature (0-12) la classe 3 è quella con maggiore varianza rispetto alle altre due classi?
- Per quale feature (0-12) le tre classi sembrano meglio separate (almeno ad occhio)?

Cliccare qui per le risposte. <!--

- Feature 2 (Ceneri)
- Feature 9 (Intensità di colore)
- Feature 6 (Flavonoidi) -->

## 4. Predizione delle etichette sul test set¶

Quanto accuratamente possiamo predire la classe (1,2,3) sulla base di una singola feature? Il codice sottostante ci consente di scoprirlo.

In [11]:

    @interact( feature=IntSlider(0,0,12) )
    def test_model(feature):
        mu, var, pi = fit_generative_model(trainx, trainy, feature)
        k = 3 # Etichette 1,2,...,k
        n_test = len(testy) # Numero di osservazioni di test
        score = np.zeros((n_test,k+1))
        for i in range(0,n_test):
            for label in range(1,k+1):
                score[i,label] = np.log(pi[label]) + \
                norm.logpdf(testx[i,feature], mu[label], np.sqrt(var[label]))
        predictions = np.argmax(score[:,1:4], axis=1) + 1
        # Conteggiamo il numero di predizioni errate
        errors = np.sum(predictions != testy)
        print("Errori di test sulla base della feature " + featurenames[feature] + ": " + str(errors) + "/" + str(n_test))

### Un ultimo esercizio¶

In questo quaderno, abbiamo considerato classificatori che usano solo una delle 13 possibili feature. La scelta di un sottoinsieme di feature è detta feature selection. In generale, questo va fatto solamente sulla base del training set (o meglio, sulla base di un apposito validation set) -- in particolare, senza considerare il test set.

Per i dati sul vino, calcolate l'errore di training e l'errore di test associato con ogni possibile scelta di una feature.

In [12]:

    ### Scrivere il proprio codice qui

Sulla base dei vostri risultati, rispondere alle seguenti domande:

- Quali tre feature hanno l'errore più basso sul training set? Elencarle dalla migliore alla peggiore.
- Quali tre feature hanno l'errore più basso sul test set? Elencarle dalla migliore alla peggiore.

Cliccare qui per la risposta. <!--

- Feature 6 (27 errori), 12 (35), 9 (38)

- Feature 6 (8 errori), 9 (10), 10 (14) Il codice utilizzato è riportato di seguito. -->

In [ ]: