# -*- coding: utf-8 -*-
"""Progetto_Data_Intensive.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1d5-w-eV5cf4S4Qf0MmpTzamtBUvwaTeK

<font size="5">Programmazione di Applicazioni Data Intensive 2020/2021 - Progetto

- Romagnoli Giacomo

- Parrinello Angelo

# Introduzione al problema

Ci troviamo sempre più spesso nella condizione di affidarci al giudizio di esperti per decretare la qualità di un vino. Fortunatamente il consumatore, scarsamente dotato in gusto e olfatto, può affidarsi a metri di giudizio scientifici per rendersi autonomo nella scelta.

Il progetto ha come **obbiettivo** quello di definire la qualità (buona o scadente) di un vino in base alle sue caratteristiche fisiche. 

In materia il problema è definito di **classificazione**.

#Caricamento Librerie

Per prima cosa carichiamo le librerie per effettuare operazioni sui dati


*   **Numpy** per creare e operare su array a N dimensioni
*   **Pandas** per caricare e manipolare dati tabulari
*   **Matplotlib** per creare grafici

Importiamo le librerie usando i loro alias convenzionali e abilitando l'inserimento dei grafici inline.
"""

# Commented out IPython magic to ensure Python compatibility.
import csv
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt

# %matplotlib inline

"""#Descrizione del problema

- Carichiamo dalla repo di Github il dataset contenente la matrice dei vini con le loro caratteristiche.
"""

import os

wine_zip = "wine.zip"

if not os.path.exists(wine_zip):
    from urllib.request import urlretrieve
    urlretrieve("https://github.com/AngeloParrinello/Wine-Quality-Classification/blob/main/wine.zip?raw=true", wine_zip)

from zipfile import ZipFile

!unzip "wine.zip"

"""- Usiamo la funzione `read_csv` per importare i dati contenuti nel file "wine.csv", specificando il separatore."""

data = pd.read_csv("wine.csv", delimiter=";")

"""- Analizziamo ora i dati, partendo dalla dimensione della matrice."""

data.shape

"""- Le colonne, ovvero le features, da cui è composta."""

data.columns

"""- Ed una piccola sezione della matrice appena importata: in questo caso visualizziamo 5 dati partendo dal fondo."""

data.tail()

"""Come si nota, il dataset contiene 6440 righe e 12 colonne.
Ogni riga costituisce un *vino registrato*. Questi dati sono stati raccolti da un'azienda vinicola portoghese. 
Di ogni vino abbiamo le seguenti *feature fisiche*:

- **fixed acidity** : *acidità fissa*. Misurato in g/L.
- **volatile acidity** : *acidità volatile*. Misurato in g/L.
- **citric acid** : *acido citrico*. Misurato in g/L.
- **residual sugar** : *zucchero residuo*. Misurato in g/L.
- **chlorides** : *cloruri*. Misurato in g/L.
- **free sulfur dioxide** : *anidride solforosa libera*. Misurata mg/L.
- **total sulfur dioxide** : *anidride solforosa totale*. Misurata mg/L.
- **density** : *densità*. Misurata g/mL.
- **pH** : *pH*.
- **sulphates** : *solfati*. Misurata in g/L.
- **alcohol** : *alcool*. Misurata in percentuale, sul totale della soluzione.

Utilizziamo il metodo *info* per vedere la memoria occupata dal DataFrame e i tipi delle colonne assegnati da Pandas, al fine di conoscere il tipo delle variabili con cui lavoreremo.
"""

data.info()

"""Per le finalità del progetto la colonna **Quality** fornisce informazioni eccessive. E' quindi necessario operare una trasformazione su quest'ultima. Nello specifico sostituiremo il valore intero attuale con un booleano: se la qualità è maggiore o uguale a 6 allora verrà etichettato come "Good" altrimenti come "Bad"; assumendo due valori, dunque, possiamo dire che è una variabile categorica e possiamo convertirla nel tipo appropriato risparmiando memoria, anche se è un dataset di piccole dimensioni. Anche la variabile *alcohol* è di tipo object e anch'essa va trasformata, questa volta in float.

Andiamo ora a mappare la variabile da predire ovvero *quality*.

- Creiamo quindi una banale funzione che adempi questo scopo.
"""

def convert_value(i):
  if i >= 6:
    return 'Good'
  else:
    return 'Bad'

"""- Attraverso la funzione *apply* di Pandas, applichiamo la funzione appena creata alla colonna **Quality**."""

data.quality = data.quality.apply(convert_value)

data['quality'].unique()

data.tail()

"""- Una volta mappata *quality*, rendiamola del corretto tipo. Così come *alcohol*."""

data = data.astype({'quality': 'category', 'alcohol': 'float64'})

"""- Tutti gli altri attributi sono variabili continue di tipo float e notiamo anche che non abbiamo valori nulli."""

data.info()

"""#Analisi Esplorativa dei Dati

Iniziamo ora *l'analisi esplorativa dei dati* ove questi vengono analizzati, esplorati e si apprendono le prime informazioni generiche a riguardo, che guideranno i passi successivi. In questa fase cercheremo anche eventuali correlazioni tra essi e risolveremo imprecisioni del dataset.

- Utilizziamo la funzione *describe*() per iniziare a comprendere come sono distribuiti i dati.

- Per evitare disambiguità al lettore precisiamo che la riga con scritto 25% (25-esimo), 50% (50-esimo) e 75% (75-esimo) stanno ad indicare il *percentile*. In breve, considerando la prima colonna e la riga del 25-esimo percentile, 6.4 è il valore sotto il quale si trovano il 25% di osservazioni.
"""

data.describe()

"""- Per analizzare la distribuzione dei valori nella colonna nominale *quality*, possiamo utilizzare *value_counts*, specificando anche *normalize=True*, se vogliamo visualizzare la percentuale."""

data["quality"].value_counts()

data["quality"].value_counts(normalize=True)

data["quality"].value_counts().plot.pie();

"""- Dai grafici possiamo notare che le classi sono **sbilanciate**, cioè sono presenti molti più vini buoni rispetto a quelli cattivi. Questa peculiarità del dataset potrebbe creare problemi in fase di learning del modello: andranno bilanciate le classi.

- Visualizziamo la distribuzione dei valori delle feature del nostro dataset in degli **istogramma**.
"""

import warnings
warnings.filterwarnings("ignore")

plt.figure(figsize=(30, 20));
for n, value in enumerate(["fixed acidity", "volatile acidity", "citric acid", "residual sugar",
                           "chlorides", "free sulfur dioxide", "total sulfur dioxide", "density", "pH", "sulphates", "alcohol"], start=1):
  data[value].plot.hist(ax=plt.subplot(4,3,n), bins=30);
  plt.subplot(4, 3, n).set_title(value);

"""- Per maggiore dettaglio gli istogrammi, ma anche altri grafici in questo notebook, saranno ripetuti per ogni feature del dataset.

- Possiamo notare, ad esempio, che la maggior parte di vini ha una densità compresa tra 0 e 20 g/mL mentre molto pochi sono estremamente densi.
"""

#percentuale dei vini estremamente densi
(len(data[data.density > 20]) / len(data.density)) * 100

"""- Possiamo anche farci un'idea più numerica della distribuzione di queste feature. Quindi eseguiamo il **binning o discretizzazione** delle colonne numeriche, utilizzando la funzione *cut*, suddividendo così i valori in fasce di uguale ampiezza per contarne poi i valori in ciascuna. """

datino = data.columns.drop("quality")
for value in datino:
  print(pd.cut(data[value],4).value_counts()) #4 fasce uguali
  print("\n")

"""- E' facile notare come molte proprietà dei nostri vini hanno divisione dei valori simili tra loro: la distribuzione non è molto omogenea.

- Visualizziamo ora le statistiche di base sui valori delle feature mediante *boxplot*.
"""

data.info()

np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning) 
plt.figure(figsize=(16, 35))
data.boxplot(ax=plt.subplot(6,2,1), column="fixed acidity", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,2), column="volatile acidity", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,3), column="citric acid", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,4), column="residual sugar", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,5), column="chlorides", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,6), column="free sulfur dioxide", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,7), column="total sulfur dioxide", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,8), column="density", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,9), column="pH", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,10), column="sulphates", by="quality", showmeans=True);
data.boxplot(ax=plt.subplot(6,2,11), column="alcohol", by="quality", showmeans=True);

"""- **ALCOHOL**: nei vini di scarsa qualità la percentuale media di alcool è più bassa rispetto ai vini di alta qualità, che però hanno una maggiore distribuzione.
- **SULPHATES**: fra buoni e cattivi vini non notiamo una grosse differenze riguardo la media, unica discriminante la maggior presenza di outliners nei vini cattivi spostati più verso il massimo. 
- **pH**: il pH nei vini di maggior pregio è leggermente più alto, quindi un vino neutro ma tendente al basico, rispetto ai vini peggiori.
- **DENSITY**: il boxplot relativo alla densità non ci è molto di aiuto: infatti la stragrande maggioranza dei vini ha una bassissima densità, come già avevamo notato.
- **TOTAL SULFUR DIOXIDE**: l'anidride solforosa totale è tendezialmente più bassa nei vini "good". Come si vede, infatti, nei vini di scarsa qualità arriviamo a toccare i 400 mg/L.
- **FREE SULFUR DIOXIDE**: l'anidride solforosa libera è simile in entrambe le categorie di vini. Notare il picco di un dato ben oltre il baffo superiore del boxplot.
- **CHLORIDES**: osservando il grafico relativo ai cloridi notiamo che ambo le tipologie di vini che stiamo analizzando hanno diversi outliers, ma in media i cloridi sono più bassi nei vini buoni.
- **RESIDUAL SUGAR**: lo zucchero residuo è simile in entrambi i tipi di vino.
- **CITRIC ACID**: l'acido citrico è simile in entrambi i tipi di vino.
- **VOLATILE ACIDITY**: l'acidità volatile è sensibilmente inferiore nel massimo per i vini buoni.
- **FIXED ACIDITY**: l'acidità fissa è simile in entrambi i tipi di vino.

- Tra le features presentate quelle che ipotiziamo siano più rilevanti sono: l'alcool, l'anidride solforosa totale e l'acidità volatile. Questo perchè notiamo sostanziali differenze nelle due categorie di vino.

- Con l'operazione di *pivoting* possiamo sdoppiare ciscuna colonna del dataframe, suddividendo i valori relativi alle classi good e bad, al fine di visualizzare meglio quanto le variabili predittive siano correlate con la classe da predire.
"""

data.pivot(columns="quality")

"""- Mostriamo ora i grafici relativi ai valori più incisivi del dataframe pivotato."""

data.pivot(columns="quality")["alcohol"].plot.hist(bins=20, stacked=True);

data.pivot(columns="quality")["total sulfur dioxide"].plot.hist(bins=20, stacked=True);

"""- In enologia, l'anidride solforosa è utilizzata sin dalle primissime fasi della produzione del vino, a partire dal mosto fino all'imbottigliamento. Nell'usare l'anidride solforosa, è opportuno sapere che una parte di questo gas si combina con alcuni componenti del mosto o del vino, mentre la restante parte resta libera, cioè non combinata. Sarà proprio la parte libera a svolgere gli importanti effetti antiossidanti e antisettici: per questo motivo è indispensabile che l'anidride solforosa si combini il meno possibile. L'anidride solforosa combinata è comunque utile, poiché nel caso in cui la frazione libera si disperda - durante le operazioni di travaso, per esempio - una piccola parte di quella combinata si libera sostituendola. (fonte: https://www.agraria.org/viticoltura-enologia/anidride-solforosa.htm)."""

data.pivot(columns="quality")["volatile acidity"].plot.hist(bins=20, stacked=True);

"""- Si definisce acidità volatile in un vino, la quantità di acido acetico presente in un vino e viene espressa in g/l di acido acetico. E’ importante determinare l’acidità volatile del vino perché la quantità di acido acetico presente è indice dello stato di sanità dell’uva, di come procede la fermentazione e poi, in generale, dello stato di conservazione del vino. (fonte: https://www.rivistadiagraria.org/articoli/anno-2009/determinazione-dellacidita-volatile-di-un-vino/)

- Possiamo visualizzare i valori medi di tutte le feature, distinti per qualità, raggruppandoli per quest'ultimo.
"""

data_by_quality = data.groupby("quality")
data_by_quality.mean()

"""- Per valutare visivamente la correlazione tra due variabili, possiamo selezionarne alcune ed evidenziare anche le due diverse classi nel *diagramma a dispersione*, differenziando i punti ad esempio per colore mediante un dizionario che associ green ai vini buoni e rosso ai vini cattivi."""

quality_color_map = {"Good": "green", "Bad": "red"}
quality_colors = data["quality"].map(quality_color_map)

plt.figure(figsize=(16, 20))
data.plot.scatter(ax=plt.subplot(4,2,1), x="residual sugar", y="chlorides", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,2), x="chlorides", y="free sulfur dioxide", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,3), x="free sulfur dioxide", y="total sulfur dioxide", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,4), x="total sulfur dioxide", y="pH", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,5), x="pH", y="alcohol", c=quality_colors);
data.plot.scatter(ax=plt.subplot(4,2,6), x="citric acid", y="residual sugar", c=quality_colors);

chlorides = data["chlorides"]
residual_sugar = data["residual sugar"]
free_sulfur_dioxide = data["free sulfur dioxide"]
total_sulfur_dioxide = data["total sulfur dioxide"]
citric_acid = data["citric acid"]

"""- Negli terzo grafico notiamo che c'è una dipendenza fra le variabili. A conferma, calcoliamo il *coefficiente di correlazione di Pearson*."""

np.mean((free_sulfur_dioxide-free_sulfur_dioxide.mean()) * (total_sulfur_dioxide-total_sulfur_dioxide.mean())) / (free_sulfur_dioxide.std() * total_sulfur_dioxide.std())

"""- Il valore del coefficiente di Pearson è vicino a 1 (correlazione diretta): al raddoppiare di una, raddoppia anche l'altra.

- Mentre le altre variabili sono scorrelate tra loro, questo è deducibile visivamente guardando i grafi. Per dimostrarlo, calcoliamo alcuni coefficenti di Pearson di queste variabili.
"""

np.mean((chlorides-chlorides.mean()) * (residual_sugar-residual_sugar.mean())) / (chlorides.std() * residual_sugar.std())

np.mean((free_sulfur_dioxide-free_sulfur_dioxide.mean()) * (chlorides-chlorides.mean())) / (free_sulfur_dioxide.std() * free_sulfur_dioxide.std())

np.mean((citric_acid-citric_acid.mean()) * (residual_sugar-residual_sugar.mean())) / (citric_acid.std() * residual_sugar.std())

"""#Normalizzazione dei dati

Procediamo ora con la normalizzazione/standardizzazione delle variabili.

In generale le variabili coinvolte in un modello di regressione
possono utilizzare scale di valori molto diverse. Questo rende difficile l'interpretazione dei modelli che andremo ad addestrare. Una soluzione più semplice è però normalizzare i dati in modo che tutte le variabili abbiano valori in un medesimo intervallo.

- Nel nostro caso abbiamo necessità di normalizzare i dati come si nota dall'esempio sotto.
"""

data[["total sulfur dioxide", "chlorides"]]

"""- Effettueremo la standardizzazione dopo aver applicato il metodo *Hold-Out* per suddividere i dati in training e validation set in preparazione all'addestramento dei modelli.

- Prima però proviamo ad addestrare un modello d'esempio senza standardizzare le feature.
"""

from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import Lasso 
from sklearn.linear_model import Perceptron
from sklearn.model_selection import train_test_split

y = data["quality"] #variabile da predire
X = data.drop(columns=['quality']) #variabili predittive

"""- Applichiamo il metodo *Hold-Out* per partizionare casualmente i dati, riservandone i 2/3 al training set e il restante al validation set. Per partizionare definiremo il parametro *random_state* (ovvero un seed per generare divisioni casuali del set) pari a 42, un numero convenzionale.

- Evitiamo per ora di utilizzare una *pipeline*, per analizzare più agevolmente il modello.
"""

X_train, X_val, y_train, y_val = \
    train_test_split(X, y, test_size=1/3, random_state=42)

"""- Addestriamo un modello **Perceptron** d'esempio."""

model = Perceptron(random_state=42) 
model.fit(X_train, y_train)
model.score(X_val, y_val)

"""- Il nostro modello, che utilizza variabili non standardizzate, produce scarsi risultati con un tasso di accuratezza che arriva a 37,1%. Ciò accade perchè le variabili presenti presentano ordini di grandezza troppo distanti.

- Proviamo adesso a **standardizzare** le feature e vedere come ciò influisca sui nostri risultati:
"""

scaler = StandardScaler()
Xn_train = scaler.fit_transform(X_train)
Xn_val = scaler.transform(X_val)

"""- Addestriamo nuovamente un modello **Perceptron**."""

model = Perceptron(random_state=42) 
model.fit(Xn_train, y_train)

"""- Una volta addestrato il modello, possiamo trovare i valori dei pesi $\mathbf{w}$ e del bias $b$ rispettivamente negli attributi `coef_[0]` e `intercept_[0]`"""

model.coef_[0]

"""- Possiamo anche capire meglio quale variabile abbia piú peso rispetto ad un'altra."""

pd.Series(model.coef_[0], index=X_train.columns)

model.intercept_[0]

"""- Valutiamo nuovamente la bontá del nostro modello:"""

model.score(Xn_val, y_val)

"""- Dopo aver normalizzato le variabili, il nostro modello produce ottimi risultati con un accuratezza del 70%. Calcoliamo ora il numero di istanze corrette e quelle errate."""

correct_class = model.predict(Xn_val) == y_val
correct_class.value_counts()

"""- Il numero di istanze che il nostro modello interpreta correttamente é di 1505 contro le 642 errate.

#Regolarizzazione dei dati ed indagine delle feature piú rilveanti

Per individuare le feature meno importanti utilizziamo la regolarizzazione **LASSO o L1** per poi andarle eventualmente ad eliminarle.
"""

model_l1 = Perceptron(penalty='l1', alpha=0.0001, random_state=42)
model_l1.fit(Xn_train, y_train)

pd.Series(model_l1.coef_[0], index=X_train.columns)

model_l1.score(Xn_val, y_val)

"""- Come si puó notare dalla cella qui sopra, applicando una regolarizzazione LASSO il modello perde precisione. Come ulteriore dimostrazione creiamo un modello Princeptron utilizzando peró solo le feature non annullate."""

X = data[["volatile acidity", "citric acid", "chlorides", "free sulfur dioxide", "density", "pH", "sulphates", "alcohol"]]

X8_train, X8_val, y8_train, y8_val = \
    train_test_split(X, y, test_size=1/3, random_state=42)

scaler = StandardScaler()
X8n_train = scaler.fit_transform(X8_train)
X8n_val = scaler.transform(X8_val)

model_8f = Perceptron(random_state=42) 
model_8f.fit(X8n_train, y8_train)

model_8f.score(X8n_val, y8_val)

"""- E' evidente quindi che per i nostri scopi andranno considerate tutte le feature. Sará comunque importante apportare una regolarizzazione a tutte le nostre feature.

Applichiamo ora la regolarizzazione **RIDGE O L2** al nostro modello d'esempio.
"""

model_l2 = Perceptron(penalty='l2', alpha=0.0001, random_state=42)
model_l2.fit(Xn_train, y_train)

pd.Series(model_l2.coef_[0], index=X_train.columns)

model_l2.score(Xn_val, y_val)

"""- Come si vede, diminuendo il peso delle feature lo score, ovvero l'accuratezza, del modello aumenta.

- E' peró risaputo che la suddivisone di istanze tra classi sbilanciate porti a molti errori di classificazione e quindi alla generazione di un modello inadeguato: il nostro dataset soffre di questo problema. Bisogna quindi sottolineare l'inaffidabilitá degli score di accuratezza ottenuti fin'ora.

- La matrice di confusione sottostante conferma quanto appena descritto: di 1379 vini buoni presenti nel validation set solo 815 sono state etichettate correttamente.
"""

from sklearn.metrics import confusion_matrix
import seaborn as sns

plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con L1', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, model_l1.predict(Xn_val)), columns=model_l1.classes_, index = model_l1.classes_), annot=True, cmap='Blues', fmt='d');

"""#Generazione di modelli di learning

##Risoluzione problema classi sbilanciate

Come dimostrato nella fase precedente c'é un forte sbilanciamento delle classi.
E' necessario quindi attuare un meccanismo di **oversampling**, ovvero aumentare il numero delle istanze nella classe meno rappresentata. Come spiegato a lezione, si utilizzerà la tecnica **SMOTE** o **Synthetic Minority Oversampling Technique** che fa parte del modulo imblearn.

- Procediamo quindi a fare l'oversampling del training set e del validation set. Come da documentazione si userà il metodo *fit_resample(X,y)*.
"""

from imblearn.over_sampling import SMOTE
import warnings
warnings.filterwarnings('ignore')

sm = SMOTE(random_state=42)

X_bal, y_bal = sm.fit_resample(X, y)
X_bal = pd.DataFrame(X_bal, columns=X.columns)
y_bal = pd.Series(y_bal, dtype="category", name="quality")

X_train_bal, y_train_bal = sm.fit_resample(X_train, y_train)
X_train_bal = pd.DataFrame(X_train_bal, columns=X_train.columns)
y_train_bal = pd.Series(y_train_bal, dtype="category", name="quality")

"""Prima del resampling:"""

pd.Series(y_train).value_counts()

pd.Series(y_train).value_counts().plot.pie(autopct='%1.1f%%',
        shadow=True, startangle=90);

"""Dopo il resampling:"""

y_train_bal.value_counts()

y_train_bal.value_counts().plot.pie(autopct='%1.1f%%',
        shadow=True, startangle=90);

"""- Le classi sono ora bilanciate.

- Creiamo prima dei dizionari per salvarci alcuni dati sui modelli che andremo a creare:
"""

scores = {}
recall = {}
precision = {}
f1 = {}
models = {}
confusion_matrices = {}

"""## Perceptron

La fase di modellazione ha inizio provando un modello basato su Perceptron. Il _Perceptron_ è un algoritmo di apprendimento molto semplice, concettualmente simile alla discesa gradiente. Per trovare i migliori parametri possibili usiamo _GridSearch_: la grid search testerá tutte le combinazioni di iperparametri possibili. Per questo modello Perceptron cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
- standardizzazione o meno delle feature
- tipo di regolarizzazione del modello
- peso della regolarizzazione
- stima o meno dell'intercetta
"""

from sklearn.model_selection import KFold
kf = KFold(5, shuffle=True, random_state=42)

per_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("per", Perceptron(n_jobs=-1, random_state=42)) #n_jobs=-1 per utilizzare tutti i processori
])

#In generale, tutti i parametri disponibili per l'ottimizzazione in GridSearchCV sono disponibili tramite
#per_mod.get_params().keys()

#ATTENZIONE! Nel dizionario per_grid (per cercare la migliore ottimizzazione 
#dei vari parametri) la chiave di ogni possibile valore del parametro dovrá essere
#esattamente la stessa contenuta in per_mod.get_params().keys()

per_grid = {
    "scaler": [None, StandardScaler()],
    "per__penalty": ["l2", "l1", "elasticnet"],
    "per__alpha": [0.0001, 0.001, 0.01, 1, 10],
    "per__fit_intercept": [False, True]
}

from sklearn.model_selection import GridSearchCV

#se voglio dentro cv_results_ i training score --> return_train_score=True 
per_gs = GridSearchCV(estimator=per_mod, param_grid=per_grid, cv=kf, n_jobs=-1)

per_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(per_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',per_gs.best_params_)
print('\n')
print('Best estimator:',per_gs.best_estimator_)

"""Oltre all'accuratezza come percentuale di classificazioni corrette, é importante tenere conto anche di altri parametri che indicano quanto il nostro modello sia preciso.

Oltre alla matrice di confusione vista prima ora utilizzeremo come metro di giudizio per la precisione del modello anche:
- F1 Score, misura unica della performance di un modello è la **_F1-measure_**, ovvero la media armonica tra precision e recall.
- Precision Score, la **_precision_** indica la percentuale di esempi classificati come "Good" che sono realmente tali.
- Recall score, la **_recall_** indica la percentuale di esempi realmente di classe "Bad" che sono stati rilevati essere tali dal modello.
"""

from sklearn.metrics import precision_score, recall_score, f1_score

#col parametro average si indica se si vuole la media(macro) o un vettore (None) o altro
#col parametro pos_label si indica la classe di riferimento

print('F1 score: {:.5f}'.format(f1_score(y_val, per_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, per_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, per_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Perceptron', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, per_gs.predict(X_val)), columns=per_gs.classes_, index = per_gs.classes_), 
            annot=True, cmap='Oranges', fmt='d');

"""Salviamo le informazioni nei vari dizionari, che più tardi ci serviranno."""

scores["per"] = per_gs.score(X_val, y_val)
recall["per"] = recall_score(y_val, per_gs.predict(X_val), pos_label='Good')
precision["per"] = precision_score(y_val, per_gs.predict(X_val), pos_label='Good')
f1["per"] = f1_score(y_val, per_gs.predict(X_val), pos_label='Good')
models["per"] = per_gs.best_estimator_
confusion_matrices["per"] = pd.DataFrame(confusion_matrix(y_val, per_gs.predict(X_val)), columns=per_gs.classes_, index = per_gs.classes_)

"""Dalla matrice di confusione e dagli score si puó notare come il modello appare ancora piuttosto sbilanciato e leggermente preciso.

Testiamo ora un modello basato sulla **Regressione Logistica**.

## Regressione Logistica

La fase di modellazione prosegue testando un modello basato sulla Regresione Logistica. La  _Regressione Logistica_ è un  un modello di classificazione binaria basato sulla regressione lineare. Per trovare i migliori parametri possibili usiamo nuovamente _GridSearch_.

Per questo modello cerchiamo i seguenti migliori iperparametri in due casi:\
*Il primo, con regolarizzazione L2 e L1 senza specificare l1_ratio:*
- standardizzazione o meno delle feature
- tipo di regolarizzazione del modello
- peso della regolarizzazione
- stima o meno dell'intercetta

*Il secondo, con elasticnet specificando i valori possibili di l1_ratio:*
- standardizzazione o meno delle feature
- tipo di regolarizzazione del modello
- peso della regolarizzazione
- peso l1_ratio
- stima o meno dell'intercetta

P.S: in Python il secondo iperparametro che la Regressione Elastic Net introduce (&alpha;) si identifica con _l1_ratio_
"""

from sklearn.linear_model import LogisticRegression
#random_state indica il seed per la casualità usato nell'addestramento
#solver indica quale implementazione utilizzare tra per l'addestramento, "saga" 
#è l'unica che supporta tutte le opzioni per la regolarizzazione

log_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("log", LogisticRegression(solver='saga', random_state=42))
])

#np.logspace come in matlab, estremi e poi numero di passi tra gli estremi
log_grid = [
    {
        "scaler": [None, StandardScaler()],
        "log__penalty": ["l2", "l1"],
        "log__C": np.logspace(-3, 3, num=10),
        "log__fit_intercept": [False, True]
    },
    {
        "scaler": [None, StandardScaler()],
        "log__penalty": ["elasticnet"],
        "log__C": np.logspace(-3, 3, num=10),
        "log__l1_ratio": np.logspace(0,1,num=10),
        "log__fit_intercept": [False, True]
    }
]

log_gs = GridSearchCV(estimator=log_mod,param_grid=log_grid, cv=kf, n_jobs=-1)

log_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(log_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',log_gs.best_params_)
print('\n')
print('Best estimator:',log_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, log_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, log_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, log_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Regressione Logistica', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, log_gs.predict(X_val)), columns=log_gs.classes_, index = log_gs.classes_), 
            annot=True, cmap='Reds', fmt='d');

scores["log"] = log_gs.score(X_val, y_val);
recall["log"] = recall_score(y_val, log_gs.predict(X_val), pos_label='Good')
precision["log"] = precision_score(y_val, log_gs.predict(X_val), pos_label='Good')
f1["log"] = f1_score(y_val, log_gs.predict(X_val), pos_label='Good')
models["log"] = log_gs.best_estimator_
confusion_matrices["log"] = pd.DataFrame(confusion_matrix(y_val, log_gs.predict(X_val)), columns=log_gs.classes_, index = log_gs.classes_)

"""Come si vede sopra il modello appena generato é piú preciso sotto tutti i punti di vista rispetto al classico Perceptron. Rimane peró leggermente alta la classificazione di falsi positivi. La modellazione prosegue con un modello di tipo **SVM** o **Support Vector Machines**.

## Support Vector Machines

Continuiamo i nostri test utilizzando ora un modello basato sul Support Vector Machines o SVM.  Il  _SVM o Support Vector Machines_ è un  un modello, definito in origine per problemi di classificazione binaria basato, con l'obiettivo di inviduare la separazione lineare ottimale tra le istanze delle due classi. Per trovare i migliori parametri possibili usiamo nuovamente _GridSearch_.

Per questo modello cerchiamo i seguenti migliori iperparametri in due casi:\
*Il primo, che utilizza la funzione Kernel Gaussiana RBF:*
- tipologia funzione Kernel
- peso delle classi
- peso della regolarizzazione
- peso del coefficente Kernel gamma (&gamma;)

*Il secondo, SVC lineare:*
- tipologia funzione Kernel
- peso delle classi
- peso della regolarizzazione

Nota: é stato tentato anche il caso con la Funzione Kernel Polinomiale ma con scarsi risultati ed eccessivi tempi di addestramento. A tal proposito non verrá riportato in questo documento.
"""

from sklearn.svm import SVC

svc_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("svc", SVC(random_state=42))
])

#svc_mod.get_params().keys()
#gamma si poteva impostare in maniera automatica con gamma='scale' quando si crea l'SVC

svc_grid =[ {
    "svc__kernel": ["rbf"],
    "svc__class_weight": [None, "balanced"],
    "svc__C": [0.01, 0.1, 1, 10,100],
    'svc__gamma': [0.01, 0.1, 1, 10,100]
}, {
    "svc__kernel": ["linear"],
    "svc__class_weight": [None, "balanced"],
    "svc__C": [0.01, 0.1, 1, 10,100],
}
]

svc_gs = GridSearchCV(estimator=svc_mod,param_grid=svc_grid, cv=kf, n_jobs=-1)

svc_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(svc_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',svc_gs.best_params_)
print('\n')
print('Best estimator:',svc_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, svc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, svc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, svc_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con SVC', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, svc_gs.predict(X_val)), columns=svc_gs.classes_, index = svc_gs.classes_), 
            annot=True, cmap='Greens', fmt='d');

scores["svc"] = svc_gs.score(X_val, y_val);
recall["svc"] = recall_score(y_val, svc_gs.predict(X_val), pos_label='Good')
precision["svc"] = precision_score(y_val, svc_gs.predict(X_val), pos_label='Good')
f1["svc"] = f1_score(y_val, svc_gs.predict(X_val), pos_label='Good')
models["svc"] = svc_gs.best_estimator_
confusion_matrices["svc"] = pd.DataFrame(confusion_matrix(y_val, svc_gs.predict(X_val)), columns=svc_gs.classes_, index = svc_gs.classes_)

"""## Albero Decisionale

Il prossimo modello che vedremo é generato attraverso un Albero Decisionale o Decision Tree. I modelli di classificazione visti finora si basano su iperpiani descritti da equazioni (lineari o non) su tutte le variabili. L' **Albero Decisionale** o **Decision Tree** costituisce un approccio differente: la classificazione avviene in base ad una serie di decisioni "semplici", basate ciascuna su una sola variabile. L'obiettivo é quindi:
- creare una struttura ad albero dove ogni nodo intermedio contiene un predicato con una variabile di input
- il predicato divide le istanze di input in 2 sottoinsiemi a cui corrispondono 2 nodi figli con relativi sotto alberi e così via ricorsivamente
"""

from sklearn.tree import DecisionTreeClassifier

features_num = X.columns.size

tree_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("tree", DecisionTreeClassifier(class_weight="balanced", random_state=42))
])

#print(tree_mod.get_params().keys())

tree_grid = {'scaler': [None, StandardScaler()],
             'tree__min_samples_split': range(2, 6),
             'tree__min_samples_leaf': range(1, 6), 
             'tree__max_depth': range(2,6),
             'tree__max_features': range(2, features_num)}

tree_gs = GridSearchCV(estimator=tree_mod,param_grid=tree_grid, cv=kf, n_jobs=-1)

tree_gs.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(tree_gs.score(X_val, y_val)))
print('\n')
print('Best parameter:',tree_gs.best_params_)
print('\n')
print('Best estimator:',tree_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, tree_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, tree_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, tree_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Albero Decisionale', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, tree_gs.predict(X_val)), columns=tree_gs.classes_, index = tree_gs.classes_), 
            annot=True, cmap='Purples', fmt='d');

"""Per visualizzare maggiormente il modello possiamo stampare _testualmente_ e _graficamente_ il tutto. Per una maggior comprensione la profonditá é stata limitata a 2."""

from sklearn.tree import export_text, plot_tree
features_name = ['fixed acidity', 'volatile acidity', 'citric acid', 'residual sugar',
       'chlorides', 'free sulfur dioxide', 'total sulfur dioxide', 'density',
       'pH', 'sulphates', 'alcohol']
print(export_text(tree_gs.best_estimator_[1], feature_names=features_name, max_depth=3))


plt.figure(figsize=(18, 12))
plot_tree(tree_gs.best_estimator_[1], feature_names=features_name,  max_depth=2, filled=True);

"""La rappresentazione testuale mostra come il modello classifichi ciascuna richiesta: 
- se l'alcohol non é maggiore di 10,79 allora considera l'aciditá volatile
  - se non superiore a 0,20 allora considera il diossido di solfuro libero altrimenti ricontrolla l'aciditá volatile.
- E cosí via fino a quando non classifica la richiesta in una maniera o in un'altra.
"""

scores["tree"] = tree_gs.score(X_val, y_val);
recall["tree"] = recall_score(y_val, tree_gs.predict(X_val), pos_label='Good')
precision["tree"] = precision_score(y_val, tree_gs.predict(X_val), pos_label='Good')
f1["tree"] = f1_score(y_val, tree_gs.predict(X_val), pos_label='Good')
models["tree"] = tree_gs.best_estimator_
confusion_matrices["tree"] = pd.DataFrame(confusion_matrix(y_val, tree_gs.predict(X_val)), columns=tree_gs.classes_, index = tree_gs.classes_)

"""## Generazione Modelli Opzionali

Per interesse personale sono stati generati alcuni modelli solamente accennati o brevemente studiati durante il corso di Data Intensive. Per tali ragioni non verrano presi in considerazione nella decisione per il miglior modello generato.

Alcune difficoltà incontrate nello sviluppo di questi modelli sono state superate informandoci anche online tra cui dalla seguente fonte https://machinelearningmastery.com/ .

### KNeighborsClassifier

*KNN* è un algoritmo di apprendimento automatico supervisionato che può essere utilizzato per risolvere problemi di classificazione e regressione. Il principio di KNN è che il valore o la classe di un dato è determinato dai dati attorno a questo valore.

L'algoritmo di previsione calcola la distanza dal punto sconosciuto x, a tutti i punti in input. I punti vengono quindi ordinati aumentando la distanza da x. La previsione viene effettuata prevedendo l'etichetta di maggioranza dai punti più vicini "K".

La scelta di una K influenzerà la classe a cui verrà assegnato un nuovo punto.

P.S: Il valore passato a knc_n_neighbors rappresenta il valore K.
"""

from sklearn.neighbors import KNeighborsClassifier

knc_mod = Pipeline([
    ('scaler', StandardScaler()),
    ('knc', KNeighborsClassifier(n_jobs=-1))
])

knc_grid = {"scaler": [None, StandardScaler()],
            'knc__n_neighbors': range(1, 10, 1),
            'knc__weights': ['uniform', 'distance']}

knc_gs = GridSearchCV(estimator=knc_mod,param_grid=knc_grid, cv=kf, n_jobs=-1)

knc_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(knc_gs.best_score_))
print('\n')
print('Best parameter:',knc_gs.best_params_)
print('\n')
print('Best estimator:',knc_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, knc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, knc_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, knc_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con KNeighborsClassifier', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, knc_gs.predict(X_val)), columns=knc_gs.classes_, index = knc_gs.classes_), 
            annot=True, cmap='flare', fmt='d');

"""### XGBoost

L’algoritmo *XGBoost*, acronimo di eXtreme Gradient Boosting, è un’implementazione specifica del metodo Gradient Boosting che utilizza approssimazioni più accurate per trovare il miglior modello ad albero. 

La formazione di un XGBoost è una procedura iterativa che calcola ad ogni passo la migliore suddivisione possibile per il k-esimo albero elencando tutte le possibili strutture ancora disponibili in quel punto del percorso. Per questo modello XGBoost cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
- xgb__eta o learning rate
- xgb__max_depth o numero massimo di nodi consentiti dalla radice alla foglia più lontana di un albero
- xgb__n_estimators o numero di alberi in un modello XGBoost
- xgb__alpha o il peso della regolarizzazione
"""

from xgboost import XGBClassifier

xgb_mod = Pipeline([
    ('scaler', StandardScaler()),
    ('xgb', XGBClassifier(n_jobs=8, random_state=42, objective='binary:logistic'))
])

xgb_grid = {
    'xgb__eta': [0.002, 0.1, 0.5],
    'xgb__max_depth': [6],
    'xgb__n_estimators': [150, 300],
    'xgb__alpha': [0.0001, 0.001]
}

xgb_gs = GridSearchCV(estimator=xgb_mod,param_grid=xgb_grid, cv=kf, n_jobs=-1)

xgb_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(xgb_gs.best_score_))
print('\n')
print('Best parameter:',xgb_gs.best_params_)
print('\n')
print('Best estimator:',xgb_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, xgb_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, xgb_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, xgb_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con XGBoost', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, xgb_gs.predict(X_val)), columns=xgb_gs.classes_, index = xgb_gs.classes_), 
            annot=True, cmap='YlOrBr', fmt='d');

"""### Random Forest

Ogni modello utilizzato dalla previsione del *Random Forest* è un albero decisionale.

Una Random Forest combina molti alberi decisionali in un unico modello. Individualmente, le previsioni fatte dai singoli alberi decisionali potrebbero non essere accurate, ma combinate insieme, le previsioni saranno in media più vicine al risultato. Per questo modello Random Forest cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
- rf__n_estimators o numero di alberi in un modello Random Forest ovvero nella foresta
- rf__max_depth o numero massimo di nodi consentiti dalla radice alla foglia più lontana di un albero
- rf__min_samples_leaf o numero minimo di campioni richiesti per trovarsi in un nodo foglia.
- rf__min_samples_split o numero minimo di campioni richiesti per dividere un nodo interno
"""

from sklearn.ensemble import RandomForestClassifier

rf_mod = Pipeline([
    ('scaler', StandardScaler()),
    ('rf', RandomForestClassifier(n_jobs=-1, random_state=42))
])

rf_grid = {
    'rf__n_estimators': [100, 200, 300],
    'rf__max_depth': [2, 4, 6, 8, 10],
    'rf__min_samples_leaf': [1, 2, 4],
    'rf__min_samples_split': [2, 5, 10],
}

rf_gs = GridSearchCV(estimator=rf_mod,param_grid=rf_grid, cv=kf, n_jobs=-1)

rf_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(rf_gs.best_score_))
print('\n')
print('Best parameter:',rf_gs.best_params_)
print('\n')
print('Best estimator:',rf_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, rf_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, rf_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, rf_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con Random Forest', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, rf_gs.predict(X_val)), columns=rf_gs.classes_, index = rf_gs.classes_), 
            annot=True, cmap='Purples', fmt='d');

"""### MLPClassifier

Lo scopo di un classificatore lineare è arrivare ad avere un iperpiano, determinato dalla combinazione lineare di $n$ variabili in input, in grado di dividere le istanze in classi. Se introducessimo delle variabili nascoste (hidden) in modo che l'iperpiano individuato dal modello sia combinazione lineare delle variabili nascoste e ciascuna nascosta sia combinazione lineare di quelle in input, allora otterremmo un *Multi Layer Perceptron* in quanto si tratta di piu modelli lineari disposti a strati. Per questo modello Multi Layer Perceptron cerchiamo di ottenere il massimo rendimento possibile dai seguenti iperparametri:
  - mlp__hidden_layer_sizes o numero di variabili nascoste da introdurre
  - mlp__batch_size o batch size (dimensione del lotto) cioè il numero di istanze che prenderà ad ogni iterazione dell'addestramento
"""

from sklearn.neural_network import MLPClassifier
MLP_mod = Pipeline([
    ("scaler", StandardScaler()),
    ("mlp", MLPClassifier(random_state=42, activation="relu"))
])

MLP_grid = {
    "mlp__hidden_layer_sizes": [16, 32, (16, 8)],
    "mlp__batch_size": [100, 200]
}

MLP_gs = GridSearchCV(estimator=MLP_mod,param_grid=MLP_grid, cv=kf, n_jobs=-1)

MLP_gs.fit(X_train_bal, y_train_bal)

print('Best CV score: {:.5f}'.format(MLP_gs.best_score_))
print('\n')
print('Best parameter:',MLP_gs.best_params_)
print('\n')
print('Best estimator:',MLP_gs.best_estimator_)

print('F1 score: {:.5f}'.format(f1_score(y_val, MLP_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, MLP_gs.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, MLP_gs.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello con MLP', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, MLP_gs.predict(X_val)), columns=MLP_gs.classes_, index = MLP_gs.classes_), 
            annot=True, cmap='Purples', fmt='d');

"""# Valutazione modelli generati

## Valutazione modelli attraverso metriche standard

Iniziamo ora a valutare i modelli fin'ora generati sulla base di alcuni informazioni estrapolate da essi.

La prima metrica che utilizzeremo per confrontare la bontà dei modelli è lo **scarto quadratico medio dell'errore dei relativi iperparametri trovati con la GridSearch**: esso corrisponde alla media dei quadrati delle differenze tra ciascun valore reale e la corrispondente predizione.
"""

from sklearn.model_selection import cross_val_score

name_model = {
    'per': 'Perceptron',
    'log': 'LogisticRegression',
    'svc': 'Support Vector Machines',
    'tree': 'Decision Tree'
}

for model in models.values():
  print(name_model[list(model.named_steps.keys())[1]] +
        ": %0.6f di accuratezza con deviazione standard di %0.6f" % 
        (cross_val_score(models[list(model.named_steps.keys())[1]],X_val, y_val,cv=kf).mean(), 
         cross_val_score(models[list(model.named_steps.keys())[1]],X_val, y_val,cv=kf).std()))
  print('\n')

"""Come possiamo notare in tutti i metodi i valori degli score sono molto vicini alla media (deviazioni standard piccole).

Confrontiamo i valori  di altri modi per valutare l'accuratezza di un classificatore: 
- **(R$^2$) Score**
- **Recall**
- **Precision**
- **F1 Score**
"""

pd.DataFrame.from_dict(scores, orient="index", columns=["(R^2) Score"])

pd.DataFrame.from_dict(recall, orient="index", columns=["Recall"])

pd.DataFrame.from_dict(precision, orient="index", columns=["Precision"])

pd.DataFrame.from_dict(f1, orient="index", columns=["F1 Score"])

"""Dai risultati emerge come tutti i modelli generati diano buoni risultati, riuscendo a classificare correttamente la maggior parte delle istanze. In particolare è chiaro come il modello basato sul SVM sia il migliore (seguito dal Decision Tree) mentre quello istanziato attraverso Perceptron sia il peggiore: ragionevole in quanto il Perceptron è l'algoritmo più semplice e vecchio tra questi testati.

A conferma di quanto appena detto osserviamo le matrici di confusione di questi:
"""

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con Perceptron', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, per_gs.predict(X_val)), columns=per_gs.classes_, index = per_gs.classes_), 
            annot=True, cmap='Oranges', fmt='d');

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con Regressione Logistica', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, log_gs.predict(X_val)), columns=log_gs.classes_, index = log_gs.classes_), 
            annot=True, cmap='Reds', fmt='d');

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con SVC', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, svc_gs.predict(X_val)), columns=svc_gs.classes_, index = svc_gs.classes_), 
            annot=True, cmap='Greens', fmt='d');

plt.figure(figsize=(4, 3))
plt.title('Matrice di confusione sul modello con Albero Decisionale', size=10)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, tree_gs.predict(X_val)), columns=tree_gs.classes_, index = tree_gs.classes_), 
            annot=True, cmap='Purples', fmt='d');

"""## Valutazione modelli con Intervallo di confidenza

Produciamo ora un confronto con intervallo di confidenza per quantificare l'accuratezza dei modelli. Fissiamo quindi la percentuale di certezza al 95%.
"""

def confidence(acc, N, Z=1.96):
    den = (2*(N+Z**2))
    var = (Z*np.sqrt(Z**2+4*N*acc-4*N*acc**2)) / den
    a = (2*N*acc+Z**2) / den
    inf = a - var
    sup = a + var
    return (inf, sup)

def calculate_accuracy(conf_matrix):
    return np.diag(conf_matrix).sum() / conf_matrix.sum().sum()

name_model = {
    'per': 'Perceptron',
    'log': 'LogisticRegression',
    'svc': 'Support Vector Machines',
    'tree': 'Decision Tree'
}

for model in models.values():
  print(name_model[list(model.named_steps.keys())[1]] +
        ": Estremo inferiore %0.6f" % 
         confidence(calculate_accuracy(confusion_matrices[list(model.named_steps.keys())[1]]), len(X_val))[0])
  
  print(name_model[list(model.named_steps.keys())[1]] +
        ": Estremo superiore %0.6f" % 
         confidence(calculate_accuracy(confusion_matrices[list(model.named_steps.keys())[1]]), len(X_val))[1])
    
  print('\n')

"""## Valutazione per confronto con Modello casuale

Creiamo ora un modello casuale, utilizzando **DummyClassifier**. Questo è uno classificatore molto semplice che mette a disposizione sklearn e che può essere utile per confrontare le sue previsioni con quelle reali.
"""

from sklearn.dummy import DummyClassifier

#strategy="uniform" garantisce uniformità randomica nella predizione
rand_mod = DummyClassifier(random_state=42,strategy="uniform")
rand_mod.fit(X_train_bal, y_train_bal)

print('CV score: {:.5f}'.format(rand_mod.score(X_val, y_val)))
print('\n')
print('F1 score: {:.5f}'.format(f1_score(y_val, rand_mod.predict(X_val), pos_label='Good')))
print('\n')
print('Precision score: {:.5f}'.format(precision_score(y_val, rand_mod.predict(X_val), pos_label='Good')))
print('\n')
print('Recall score: {:.5f}'.format(recall_score(y_val, rand_mod.predict(X_val), pos_label='Good')))
print('\n')
plt.figure(figsize=(8, 6))
plt.title('Matrice di confusione sul modello Casuale', size=16)
sns.heatmap(pd.DataFrame(confusion_matrix(y_val, rand_mod.predict(X_val)), columns=rand_mod.classes_, index = rand_mod.classes_), 
            annot=True, cmap='Pastel1', fmt='d');

scores["rand"] = rand_mod.score(X_val, y_val);
recall["rand"] = recall_score(y_val, rand_mod.predict(X_val), pos_label='Good')
precision["rand"] = precision_score(y_val, rand_mod.predict(X_val), pos_label='Good')
f1["rand"] = f1_score(y_val, rand_mod.predict(X_val), pos_label='Good')
#models["rand"] = rand_mod
confusion_matrices["rand"] = pd.DataFrame(confusion_matrix(y_val, rand_mod.predict(X_val)), columns=rand_mod.classes_, index = rand_mod.classes_)

from scipy import stats

def difference_between_two_models(error1, error2, confidence):
    z_half_alfa = stats.norm.ppf(confidence)
    variance = (((1 - error1) * error1) / len(y_val)) + (((1 - error2) * error2) / len(y_val))
    d_minus = abs(error1 - error2) - z_half_alfa * (pow(variance, 0.5))
    d_plus = abs(error1 - error2) + z_half_alfa * (pow(variance, 0.5))
    print("Valore minimo: {}\nValore massimo: {}\n".format(d_minus, d_plus))

for model in models.values():
  print(name_model[list(model.named_steps.keys())[1]] + " VS Modello Random")
  print("\n")
  difference_between_two_models(1 - f1[list(model.named_steps.keys())[1]], 1 - f1["rand"], 0.99)

"""Ora abbiamo la certezza al 99% che i diversi modelli addestrati siano migliori di un modello casuale: nessuno dei valori che possiamo visualizzare contiene lo 0.

Risulta ormai evidente come il modello migliore sia quello addestrato con **Support Vector Machines**.

# Confronto finale e Conclusioni

Individuato il modello SVM come quello migliore analizziamolo.

Per individuare l'iperpiano di separazione migliore tra le due classi, SVM si basa sulle istanze più vicine ad esso, i cosiddetti **Support Vector**. Vediamo quindi quali features sono più positivamente/negativamente correlate e in che misura con la variabile *quality*.

Vediamo prima quante istanze son stato necessarie per individuare l'iperpiano:
"""

len(models['svc'][1].support_)

"""Come notiamo le istanze utilizzate son molteplici (come si nota sotto, circa il 64%), questo indica come una divisione netta e precisa tra le due classi sia difficile da individuare. Ancora più complicato per modelli lineare come Perceptron: un altro motivo per la sua scarsa precisione."""

len(models['svc'][1].support_)  / len(X_train_bal) * 100

"""Per visualizzare meglio il tutto possiamo definire una funzione che ci permetta di capire quanti punti siano stati usati per definire l'iperpiano. Plottando i punti che abbiamo usato per il training e quelli usati dal support vector, vedremo come molti combacino"""

def plot_data_and_support(X, model):

    #in Xs metto tutte le righe e colonne di X che hanno indice preso da model.support_
    Xs = X.iloc[model.support_]
    #print(Xs)

    plt.scatter(X.iloc[:, 3], X.iloc[:, 6], s=4, c='black')
    plt.scatter(Xs.iloc[:, 3], Xs.iloc[:, 6], s=4, c='yellow')

plot_data_and_support(X_train_bal, models['svc'][1])

"""Infine scopriamo quale feature son state più o meno rilevanti per i nostri scopi di classificazione."""

support_dataframe = pd.DataFrame(models['svc'][1].support_vectors_, columns=X_train_bal.columns)
#support_dataframe.describe()

mean = pd.Series([0.061903,0.035754,0.043108,-0.031508,0.083313,-0.021621,-0.036785,0.026544,0.031376,0.069159,0.016729], index= X_train_bal.columns)
mean.plot.bar()
plt.grid(axis="y")

"""Come si era previsto nella fase esplorativa dei dati, tutte le features hanno un peso più o meno significativo sul risultato del modello, è quindi importante considerarle tutte per determinare se un vino è buono o cattivo.

In definitiva considerando i valori ottenuti dal modello, ci riteniamo molto soddisfatti di questa nostra prima esperienza nel mondo del Data Science. Le valutazioni affrontate nel finale del documento ci mostrano come il modello SVM possa essere un ottimo strumento di supporto per la classificazione di un vino buono o cattivo.

E' stato scelto questo dataset dal momento che la famiglia di un membro del team di questo progetto (Romagnoli Giacomo) è proprietaria di un'azienda vinicola. Ci immaginiamo ad esempio che un'azienda come quella possa utilizzare questo dispositivo per scremare in un primo momento i vini buoni da quelli cattivi per poi farli assaggiare ad un ipotetico somelier che ne giudicherà l'effettiva bontà.
"""