# TEN SKRYPT MA BYĆ URUCHAMIANY W JUPYTER NOTEBOOK (UŻYŁEM VS CODE)

# %%
import pandas as pd
import numpy as np
from sklearn.preprocessing import PowerTransformer
from sklearn.covariance import MinCovDet
from scipy.stats import chi2

import seaborn as sb
import matplotlib.pyplot as plt



# %%
def add_is_outlier_IQR(data, col_name):
    col_values = data[col_name]
    
    Q1=col_values.quantile(0.25)
    Q3=col_values.quantile(0.75)
    IQR=Q3-Q1
    
    outliers_col_name = f'is_{col_name.replace(" ", "_")}_outlier'
    data[outliers_col_name] = ((col_values < (Q1 - 1.5 * IQR)) | (col_values > (Q3 + 1.5 * IQR)))
    
    return data


def boxPlot(data, varx, vary, title, xlab, ylab, hue = None):
    hplot = sb.boxplot(varx, vary, hue=hue, data=data)
    plt.title(title, fontsize=18)
    plt.xlabel(xlab, fontsize=16)
    plt.ylabel(ylab, fontsize=16)
    
    return hplot


def yeo_johnson_transf(data):
    pt = PowerTransformer(method='yeo-johnson', standardize=True)
    pt.fit(data)

    lambdas = pt.lambdas_
    df_yeojohnson = pd.DataFrame( pt.transform(data), columns=data.columns.values )
    
    return df_yeojohnson, lambdas


# %%
# Załaduj dane o czerwonym winie
df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv',
                 sep=';')

# Pobierz nazwy liczbowych kolumn poza kolumną opisującą jakość 
numeric_col_names = df.drop('quality', axis=1).columns.values

df

# %%
sb.set(style="whitegrid")

# Wykres zawartości siarczynów w celu sprawdzenia, czy 
# istnieją jednoznaczne wartości odstające 
boxPlot(df, varx='sulphates', vary=None,
        title='Rozkład zawartości siarczynów',
        xlab='siarczyny', ylab=None)
# W przypadku, gdy nie korzystasz z notebooka Jupyter, uruchom również następującą instrukcję: 
# plt.show()

# %%
# Ponieważ są wartości odstające, dodajmy kolumnę boolean  
# do ramki danych oznaczającą, który wiersz
# zawiera odstającą wartość odpowiadającą zawartości siarczynów
add_is_outlier_IQR(df, 'sulphates')

# Wykres pudełkowy po usunięciu początkowych wartości odstających
df_no_outliers = df.loc[~df['is_sulphates_outlier']]

boxPlot(df_no_outliers, varx='sulphates', vary=None,
        title='Rozkład zawartości siarczynów bez wartości odstających',
        xlab='siarczyny', ylab=None)
# W przypadku, gdy nie korzystasz z notebooka Jupyter, uruchom również następującą instrukcję: 
# plt.show()

# %%
# Narysujmy wykresy pudełkowe dla każdego głosowania jakości, 
# po usunięciu wstępnych wartości odstających
boxPlot(df_no_outliers, varx='quality', vary='sulphates',
        title='Rozkład zawartości siarczynów bez wartości odstających według jakości',
        xlab='jakość', ylab='siarczyny')
# W przypadku, gdy nie korzystasz z notebooka Jupyter, uruchom również następującą instrukcję: 
# plt.show()

# %%
# ANALIZA WIELOWYMIAROWA 
#-----------------------
# Wykreślenie histogramu dla wszystkich zmiennych 
# z wykorzystaniem zbioru danych bez wartości odstających
df_no_outliers.drop('quality', axis=1).hist(figsize=(10,10))
plt.tight_layout()
plt.show()

# %%
# Zastosowanie przekształceń Yeo-Johnsona w celu
# wyelimowania skośności
df_transf, lambda_arr = yeo_johnson_transf(df_no_outliers[numeric_col_names])

# Wykreślenie histogramu dla wszystkich przekształconych zmiennych 
# w celu sprawdzenia, czy skośność uległa zmniejszeniu 
df_transf.hist(figsize=(10,10))
plt.tight_layout()
plt.show()

# %%
# # OSTRZEŻENIE:  Utworzenie poniższych wykresów zajmuje kilka minut
# #
# # Jeśli chcesz również sprawdzić wykresy gęstości każdej zmiennej oraz 
# # wykresy rozproszone między poszczególnymi parami, pogrupowanymi według jakości, 
# # możesz skorzystać z wykresu par. Poniższy korzysta z wyjściowej ramki danych ... 
# sb.pairplot(df, hue='quality', diag_kind = 'kde',
#             plot_kws = {'alpha': 0.6, 's': 80, 'edgecolor': 'k'})

# # %%
# # ... natomiast ten jest wygenerowany przy użyciu przekształconej ramki danych. 
# df_transf_qual = df_transf.copy()
# df_transf_qual['quality'] = df['quality']

# sb.pairplot(df_transf_qual, hue='quality', diag_kind = 'kde',
#             plot_kws = {'alpha': 0.6, 's': 80, 'edgecolor': 'k'})
# # W przypadku, gdy nie korzystasz z notebooka Jupyter, uruchom również następującą instrukcję: 
# # plt.show()
# %%
# Obliczmy kwadratowe odległości Mahalanobisa za pomocą 
# minimalnego wyznacznika kowariancji w celu obliczenia 
# rozbudowanej macierzy kowariancji 
robust_cov = MinCovDet(support_fraction=0.7).fit(df_transf)
center = robust_cov.location_

D = robust_cov.mahalanobis(df_transf - center)
D

# %%
# Odległość kwadratów Mahalanobisa (D) jest zgodna z rozkładem chi-kwadrat 
# (https://markusthill.github.io/mahalanbis-chi-squared/#the-squared-mahalanobis-distance-follows-a-chi-square-distribution-more-formal-derivation)
#
# Biorąc pod uwagę wartość odcięcia związaną z istotnością statystyczną 
# z którą chcemy określić wartości odstające, otrzymujemy odpowiednią 
# wartość progową, powyżej której, należy uznać obserwację za wartość odstającą
cutoff = 0.98
degrees_of_freedom = df_transf.shape[1]  # przy określonej liczbie zmiennych (kolumny) 
cut = chi2.ppf(cutoff, degrees_of_freedom) # wartość progowa

# Odległość kwadratów Mahalanobisa dla wartości odstających 
D[D > cut]

# %%
# Obliczenie prawdopodobieństwa, że odległość D[5]
# jest wartością odstającą
chi2.cdf(D[5], degrees_of_freedom)

# %%
# Obliczenie, czy obserwacja jest wartością odstającą dla danego odcięcia
is_outlier_arr = (D > cut)

# Obliczenie prawdopodobieństwa, że obserwacja jest wartością odstającą, a nie przypadkiem
outliers_stat_proba = np.zeros(len(is_outlier_arr))

for i in range(len(is_outlier_arr)):
    outliers_stat_proba[i] = chi2.cdf(D[i], degrees_of_freedom)

# Ile wartości odstających o znaczeniu statystycznym większym niż granica
len(outliers_stat_proba[outliers_stat_proba > cutoff])


# %%
# Dodawanie informacji o wartości odstających do ramki danych według 
# odległości kwadratów Mahalanobisa
df['is_mahalanobis_outlier'] = is_outlier_arr
df['mahalanobis_outlier_proba'] = outliers_stat_proba

df[df['is_mahalanobis_outlier']]


# %%
