# Classificazione di vini attraverso gaussiane multivariate¶

In questo quaderno riprendiamo la classificazione dei vini, utilizzando l'intero insieme delle 13 feature.

## 1. Caricamento del dataset¶

Come al solito, iniziamo caricando il dataset.

Ricordiamo che vi sono 178 osservazioni, ognuna con 13 feature ed una etichetta (1, 2, 3). Come prima, le divideremo in un training set di 130 osservazioni ed un test set di 48 osservazioni.

In [1]:

    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    # Modulo utile per gestire distribuzioni gaussiane
    from scipy.stats import norm, multivariate_normal

In [2]:

    # Carichiamo il dataset
    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']
    # Suddividiamo le 178 istanze in un training set (trainx, trainy) di taglia 130 e un test set (testx, testy) di taglia 48
    np.random.seed(0)
    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]

## 2. Fit di un modello generativo gaussiano¶

Definiamo ora una funzione che fa il fit di un modello generativo gaussiano. Per ogni classe (`j=1,2,3`), abbiamo:

- `pi[j]`: il peso relativo della classe
- `mu[j,:]`: la media, un vettore 13-dimensionale
- `sigma[j,:,:]`: la matrice di covarianza 13x13

Ciò significa che `pi` sarà un array 4 $\times$ 1 (gli array Python sono numerati dall'indice 0, e noi lasceremo `j=0` inutilizzato), `mu` sarà un array 4 $\times$ 13 e `sigma` sarà un array 4 $\times$ 13 $\times$ 13.

In [3]:

    def fit_generative_model(x,y):
        k = 3  # numero di classi (etichette 1,2,...,k)
        d = (x.shape)[1]  # numero di feature
        mu = np.zeros((k+1,d))
        sigma = np.zeros((k+1,d,d))
        pi = np.zeros(k+1)
        for j in range(1,k+1):
            indices = (y == j) # vettore rappresentante il sottoinsieme degli esempi di classe j
            mu[j] = np.mean(x[indices,:], axis=0) # axis=0 indica che la media è ottenuta aggregando le righe,
                                                        # ovvero colonna per colonna
            sigma[j] = np.cov(x[indices,:], rowvar=0, bias=1) # rowvar=0 indica che la matrice x contiene
                                                    # le osservazioni sulle righe e le feature sulle colonne;
                                                    # bias=1 corrisponde alla normalizzazione vista a lezione
            pi[j] = float(sum(indices))/float(len(y)) # = frazione degli esempi di tipo label su tutti gli esempi
        return mu, sigma, pi

In [4]:

    # Fit di un modello generativo gaussiano sui dati di training
    mu, sigma, pi = fit_generative_model(trainx, trainy)

## 3. Predizione delle etichette per gli esempi di test¶

Esercizio da svolgere</font>: Vogliamo definire una routine di testing generale che riceva in input:

- gli array `pi`, `mu`, `sigma` che definiscono il modello generativo, come costruiti nella precedente sezione
- un test set (una matrice di input `tx` e le corrispondenti etichette `ty`)
- una lista di feature `features` (scelte nell'intervallo 0-12)

La routine deve restituire il numero di misclassificazioni del modello generativo sui dati di test, quando il modello è costruito sulla base delle sole feature specificate. Per esempio, usando solo le tre feature 2 (`'Ceneri'`), 4 (`'Magnesio'`) e 6 (`'Flavonoidi'`) si ottengono 7 misclassificazioni (su 48 esempi di test), quindi

        test_model(mu, sigma, pi, [2,4,6], testx, testy)

dovrebbe stampare `7/48`.

Per restringere l'attenzione ad un sottoinsieme di feature, scegliamo le corrispondenti coordinate del vettore 13-dimensionale delle medie e una sottomatrice opportuna della matrice di covarianza 13 $\times$ 13.

La seguente implementazione manca solo dell'istruzione che calcola il discriminante $\delta_j(x)$. Completare l'implementazione.

In [8]:

    # Testa la performance di un modello gaussiano costruito sulla base di un sottoinsieme delle feature
    def test_model(mu, sigma, pi, features, tx, ty):
        k = 3 # 3 classi
        mistakes = 0
        for i in range(len(tx)): # ciclo sugli esempi del test set
            x = tx[i] # input i-esimo
            y = ty[i] # output i-esimo
            # Ricerca del massimo discriminante: argmax_j delta_j(x)
            delta_best = -np.inf # valore massimo del discriminante trovato sinora (inizialmente è -INFINITO)
            best = None # indice corrispondente
            for j in range(1,k+1):
                mu_j = mu[j, features]
                rows = [f_row for f_row in features for f_col in features]
                cols = [f_col for f_row in features for f_col in features]
                sigma_j = sigma[j, rows, cols].reshape((len(features), len(features)))
                ### Inserire qui l'istruzione appropriata per il calcolo del discriminante:
                # delta_j = ...
                delta_j = np.log(pi[j]) + multivariate_normal.logpdf(x[features], mean=mu_j, cov=sigma_j)
                ###
                # Aggiornamento del discriminante massimo
                if delta_j > delta_best:
                    delta_best = delta_j
                    best = j
            if y != best: mistakes += 1
        print(mistakes, '/', len(tx))

Cliccare qui per la soluzione.

In [9]:

    test_model(mu, sigma, pi, [2,4,6], testx, testy)

    7 / 48

### Esercizi rapidi¶

Esercizio 1. Quante misclassificazioni si hanno sul test set utilizzando solamente la feature 2 (Ceneri)?

Cliccare qui per la risposta.

Esercizio 2. Quante misclassificazioni si hanno utilizzando le feature 0 (Alcool) e 2 (Ceneri)?

Cliccare qui per la risposta.

Esercizio 3. Quante misclassificazioni si hanno utilizzando le feature 0 (Alcool), 2 (Ceneri), e 6 (Flavonoidi)?

Cliccare qui per la risposta.

Esercizio 4. Quante misclassificazioni si hanno utilizzando tutte e 13 le feature?

Cliccare qui per la risposta.

Esercizio 5. A lezione, abbiamo ottenuto delle risposte lievemente diverse ad alcune di queste domande. Quale potrebbe essere il motivo?

Cliccare qui per la risposta.

In [ ]: