library(dplyr)
library(stringr)
library(recipes)
library(robust)


add_is_outlier_IQR <- function(data, col_name) {
    
    x <- data[[col_name]]
    quar <- quantile(x, probs = c(0.25, 0.75), na.rm = TRUE)
    iqr <- diff(quar)
    k <- 1.5
    
    outliers_col_name <- str_glue('is_{str_replace(col_name, " ", "_")}_outlier')
     
    data[[outliers_col_name]] <- ifelse((x < quar[1] - k * iqr) | (x > quar[2] + k * iqr), 1, 0)
    
    return(data)
}


yeo_johnson_transf <- function(data) {
    
    rec <- recipe(data, quality ~ .)
    
    rec <- rec %>%
        step_center( all_numeric(), - all_outcomes() ) %>%
        step_scale( all_numeric(), - all_outcomes() ) %>%
        step_YeoJohnson( all_numeric(), -all_outcomes() )
    
    prep_rec <- prep( rec, training = data )
    
    res_list <- list( df_yeojohnson = bake( prep_rec, data ),
                      lambdas = prep_rec$steps[[3]][["lambdas"]] )
}



# Pobierz nazwy liczbowych kolumn poza kolumną opisującą jakość 
numeric_col_names <- dataset %>% 
    select( where(is.numeric), -quality ) %>% 
    names()

# 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
df <- add_is_outlier_IQR(dataset, col_name = 'sulphates')


# Zastosowanie przekształceń Yeo-Johnsona w celu
# wyeliminowania skośności
yeo_johnson_list <- df %>% 
    yeo_johnson_transf()

df_transf <- yeo_johnson_list$df_yeojohnson


# Obliczmy kwadratowe odległości Mahalanobisa za pomocą 
# minimalnego wyznacznika kowariancji w celu obliczenia 
# rozbudowanej macierzy kowariancji 
data <- df_transf %>%
    select( all_of(numeric_col_names) )

cov_obj <- data %>% 
    covRob( estim="mcd", alpha=0.7 )

center <- cov_obj$center
cov <- cov_obj$cov

distances <- data %>%
    mahalanobis( center=center, cov=cov )


# Na podstawie wartości odcięcia związanej 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 <- ncol(data) # given by the number of variables (columns)

outliers_value_cutoff <- qchisq(cutoff, degrees_of_freedom) # threshold value


data <- data %>% 
    mutate(
        
        # Kolumna wskaźnika wartości odstających siarczynów
        is_sulphates_outlier      = df$is_sulphates_outlier,
        
        # Kolumna wskaźnika wartości odstających wykrytych za pomocą odległości Mahalanobisa
        is_mahalanobis_outlier    = distances > outliers_value_cutoff,
        
        # Obliczenie prawdopodobieństwa, że obserwacja jest wartością odstającą, a nie przypadkiem
        mahalanobis_outlier_proba = pchisq(distances, ncol(data))
    )

