---
title: "Penalized Discriminant Analysis on Simulated Examples"
output: html_document
---

```{r packages}
source("fun.R")
library(classPP)
library(tidyr)
library(DiscriMiner)
library(ggplot2)
library(latex2exp)
library(caret)
library(e1071)
nrun = 20
```

## Simulation 1: percentage of important variables

```{r}
lda_acc  = c()
pda4_acc = c()
pda6_acc = c()
pda8_acc = c()
svmlinear_acc = c()
svmrbf_acc = c()

x= seq(0.03, 0.05, by = 0.001) #percentage of importance variables
for (i in x){ 
  print(i)
  n = 40
  p = 500
  num_unimp = round((1-i)*p, digits = 0)
  num_imp = p - num_unimp
  
  lda_temp  = c()
  pda4_temp = c()
  pda6_temp = c()
  pda8_temp = c()
  svmlinear_temp = c()
  svmrbf_temp = c()
  
  for (j in 1:nrun){ #compute accuracy score for 20 times
    df = scale(as.data.frame(cbind(
      rbind(matrix(rnorm(n/2*num_imp,2.2),,num_imp),
            matrix(rnorm(n/2*num_imp, -2.2),,num_imp)),
      matrix(rnorm(n*num_unimp),,num_unimp))))
    class = c(rep(1, n/2), rep(2, n/2))

    ind = sample(1:40, 10, replace=FALSE)
    test = df[ind,]
    train = df[-ind,]
    cls_test = class[ind]
    cls_train = class[-ind]
  
    PP.opt = PP.optimize.anneal("LDA",2,train,cls_train)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    lda_temp = c(lda_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                    cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.4)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda4_temp = c(pda4_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.6)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda6_temp = c(pda6_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.8)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda8_temp = c(pda8_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    
    cls_train = as.factor(cls_train)
    train2 = cbind(as.data.frame(train), cls_train)
    train_ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmLinear", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmlinear_temp = c(svmlinear_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))

    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmRadial", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmrbf_temp = c(svmrbf_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))
  }
  
  lda_acc  = c(lda_acc,  mean(lda_temp))
  pda4_acc = c(pda4_acc, mean(pda4_temp))
  pda6_acc = c(pda6_acc, mean(pda6_temp))
  pda8_acc = c(pda8_acc, mean(pda8_temp))
  svmlinear_acc = c(svmlinear_acc, mean(svmlinear_temp))
  svmrbf_acc = c(svmrbf_acc, mean(svmrbf_temp))
}

write.csv(cbind(x, lda_acc, pda4_acc, pda6_acc, pda8_acc, svmlinear_acc, svmrbf_acc),
          "1_perc_imp_var.csv")
```

## Simulation 2: Ratio between dimension and observation

```{r}
lda_acc  = c()
pda4_acc = c()
pda6_acc = c()
pda8_acc = c()
svmlinear_acc = c()
svmrbf_acc = c()
x = seq(10, 100, by = 6)
for (i in x){ 
  print(i)
  n = i
  p = 500
  num_imp = 10
  num_unimp = p-num_imp
  
  lda_temp  = c()
  pda4_temp = c()
  pda6_temp = c()
  pda8_temp = c()
  svmlinear_temp = c()
  svmrbf_temp = c()
  
  for (j in 1:nrun){ #compute accuracy score for 20 times
    df = scale(as.data.frame(cbind(
      rbind(matrix(rnorm(n/2*num_imp, 2.2),,num_imp),
            matrix(rnorm(n/2*num_imp, -2.2),,num_imp)),
      matrix(rnorm(n*num_unimp),,num_unimp))))
    class = c(rep(1, n/2), rep(2, n/2))

    ind = sample(1:n, round(n*0.25, 0), replace=FALSE)
    test = df[ind,]
    train = df[-ind,]
    cls_test = class[ind]
    cls_train = class[-ind]
  
    PP.opt = PP.optimize.anneal("LDA",2,train,cls_train)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    lda_temp = c(lda_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                    cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.4)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda4_temp = c(pda4_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.6)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda6_temp = c(pda6_temp, acc(PP.opt$proj.best, t(train), t(test),
                                      cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.8)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda8_temp = c(pda8_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    
    cls_train = as.factor(cls_train)
    train2 = cbind(as.data.frame(train), cls_train)
    train_ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmLinear", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmlinear_temp = c(svmlinear_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))

    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmRadial", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmrbf_temp = c(svmrbf_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))
  }
  lda_acc  = c(lda_acc,  mean(lda_temp))
  pda4_acc = c(pda4_acc, mean(pda4_temp))
  pda6_acc = c(pda6_acc, mean(pda6_temp))
  pda8_acc = c(pda8_acc, mean(pda8_temp))
  svmlinear_acc = c(svmlinear_acc, mean(svmlinear_temp))
  svmrbf_acc = c(svmrbf_acc, mean(svmrbf_temp))
}

write.csv(cbind(x, lda_acc, pda4_acc, pda6_acc, pda8_acc, svmlinear_acc, svmrbf_acc), "2_ratio_dim_obs.csv")
```

## Simulation 3: effect of outliers

```{r}
lda_acc  = c()
pda4_acc = c()
pda6_acc = c()
pda8_acc = c()
svmlinear_acc = c()
svmrbf_acc = c()
x = seq(0, 16, by = 1)
for (i in x){ 
  n = 40
  p = 500
  num_imp = 10
  num_unimp = p-num_imp
  
    
  lda_temp  = c()
  pda4_temp = c()
  pda6_temp = c()
  pda8_temp = c()
  svmlinear_temp = c()
  svmrbf_temp = c()
  for (j in 1:nrun){ 
    df = scale(as.data.frame(cbind(
      rbind(matrix(rnorm(n/2*num_imp, 2.2),,num_imp), 
            matrix(rnorm(n/2*num_imp, -2.2),,num_imp)),
      matrix(rnorm(n*num_unimp),,num_unimp))))
  
    if (i != 0){
      #for (var in 1:num_imp){ #add outliers to all important variables
      for (var in (num_imp+1):p){ #add outliers to all variables except important variables
        for (num_outlier in 1:i){
          ind = sample.int(n,1)
          df[ind, var] = rnorm(1, 0, 15)
        }
      }
    }
  
    df = scale(df)
    class = c(rep(1, n/2), rep(2, n/2))

    ind = sample(1:n, round(n*0.25, 0), replace=FALSE)
    test = df[ind,]
    train = df[-ind,]
    cls_test = class[ind]
    cls_train = class[-ind]

    PP.opt = PP.optimize.anneal("LDA",2,train,cls_train)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    lda_temp = c(lda_temp, acc(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.4)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda4_temp = c(pda4_temp, acc(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.6)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda6_temp = c(pda6_temp, acc(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.8)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda8_temp = c(pda8_temp, acc(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    cls_train = as.factor(cls_train)
    train2 = cbind(as.data.frame(train), cls_train)
    train_ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmLinear", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmlinear_temp = c(svmlinear_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))

    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmRadial", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmrbf_temp = c(svmrbf_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))
  }
  lda_acc  = c(lda_acc,  mean(lda_temp))
  pda4_acc = c(pda4_acc, mean(pda4_temp))
  pda6_acc = c(pda6_acc, mean(pda6_temp))
  pda8_acc = c(pda8_acc, mean(pda8_temp))
  svmlinear_acc = c(svmlinear_acc, mean(svmlinear_temp))
  svmrbf_acc = c(svmrbf_acc, mean(svmrbf_temp))
}


write.csv(cbind(x, lda_acc, pda4_acc, pda6_acc, pda8_acc, svmlinear_acc, svmrbf_acc), "4_outliers.csv")
```

## Simulation 4: number of classes

```{r}
lda_acc  = c()
pda4_acc = c()
pda6_acc = c()
pda8_acc = c()
lda_mse  = c()
pda4_mse = c()
pda6_mse = c()
pda8_mse = c()
lda_mae  = c()
pda4_mae = c()
pda6_mae = c()
pda8_mae = c()
svmlinear_acc = c()
svmrbf_acc = c()
svmlinear_mse = c()
svmrbf_mse = c()
svmlinear_mae = c()
svmrbf_mae = c()

x = seq(2, 12, by = 1)
for (i in x){
  print(i)
  n = 20*i
  p = 500
  num_imp = 10
  num_unimp = p-num_imp
  
  lda_temp  = c()
  lda_mse_temp = c()
  lda_mae_temp = c()
  pda4_temp = c()
  pda4_mse_temp = c()
  pda4_mae_temp = c()
  pda6_temp = c()
  pda6_mse_temp = c()
  pda6_mae_temp = c()
  pda8_temp = c()
  pda8_mse_temp = c()
  pda8_mae_temp = c()
  svmlinear_temp = c()
  svmrbf_temp = c()
  svmlinear_mse_temp = c()
  svmrbf_mse_temp = c()
  svmlinear_mae_temp = c()
  svmrbf_mae_temp = c()
  for (j in 1:nrun){
    df = rbind(matrix(rnorm(20*num_imp, 0),,num_imp),
               matrix(rnorm(20*num_imp, 4),,num_imp))
    class = c(rep(1, n/i), rep(2, n/i))
    mean = 4
    cls = 2
    if (i!=2){
      for (j in 1:(i-2)){
        mean = 4 + mean
        cls = 1 + cls
        df = rbind(df, matrix(rnorm(20*num_imp, mean=mean),,num_imp))
        class = c(class, rep(cls, n/i))
      }
    }
    
    df = scale(as.data.frame(cbind(df, 
                                   matrix(rnorm(n*num_unimp),,num_unimp))))
  
    ind = sample(1:n, round(n*0.25, 0), replace=FALSE)
    test = df[ind,]
    train = df[-ind,]
    cls_test = class[ind]
    cls_train = class[-ind]
  
    PP.opt = PP.optimize.anneal("LDA",2,train,cls_train)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    lda_temp = c(lda_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                    cls_train, cls_test))
    lda_mse_temp = c(lda_mse_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[1])
    lda_mae_temp = c(lda_mae_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[2])
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.4)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda4_temp = c(pda4_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    pda4_mse_temp = c(pda4_mse_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[1])
    pda4_mae_temp = c(pda4_mae_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[2])
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.6)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda6_temp = c(pda6_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    pda6_mse_temp = c(pda6_mse_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[1])
    pda6_mae_temp = c(pda6_mae_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[2])
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.8)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda8_temp = c(pda8_temp, acc(PP.opt$proj.best, t(train), t(test), 
                                      cls_train, cls_test))
    pda8_mse_temp = c(pda8_mse_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[1])
    pda8_mae_temp = c(pda8_mae_temp, metrics(PP.opt$proj.best, t(train), t(test), cls_train, cls_test)[2])
    
    cls_train = as.factor(cls_train)
    train2 = cbind(as.data.frame(train), cls_train)
    train_ctrl <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmLinear", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmlinear_temp = c(svmlinear_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))
    svmlinear_mse_temp = c(svmlinear_mse_temp, mean((as.numeric(y_pred) - cls_test)^2))
    svmlinear_mae_temp = c(svmlinear_mae_temp, mean(abs(as.numeric(y_pred) - cls_test)))

    svm_classifier <- train(cls_train ~., data = train2, 
                               method = "svmRadial", 
                               trControl=train_ctrl,
                               tuneLength = 10)
    svmrbf_temp = c(svmrbf_temp, 
                    sum(predict(svm_classifier, test) == cls_test)/nrow(test))
    svmrbf_mse_temp = c(svmrbf_mse_temp, mean((as.numeric(y_pred) - cls_test)^2))
    svmrbf_mae_temp = c(svmrbf_mae_temp, mean(abs(as.numeric(y_pred) - cls_test)))
  }
  lda_acc  = c(lda_acc,  mean(lda_temp))
  pda4_acc = c(pda4_acc, mean(pda4_temp))
  pda6_acc = c(pda6_acc, mean(pda6_temp))
  pda8_acc = c(pda8_acc, mean(pda8_temp))
  lda_mse  = c(lda_mse,  mean(lda_mse_temp))
  pda4_mse = c(pda4_mse, mean(pda4_mse_temp))
  pda6_mse = c(pda6_mse, mean(pda6_mse_temp))
  pda8_mse = c(pda8_mse, mean(pda8_mse_temp))
  lda_mae  = c(lda_mae,  mean(lda_mae_temp))
  pda4_mae = c(pda4_mae, mean(pda4_mae_temp))
  pda6_mae = c(pda6_mae, mean(pda6_mae_temp))
  pda8_mae = c(pda8_mae, mean(pda8_mae_temp))
  svmlinear_mse = c(svmlinear_mse, mean(svmlinear_mse_temp))
  svmrbf_mse = c(svmrbf_mse, mean(svmrbf_mse_temp))
  svmlinear_mae = c(svmlinear_mae, mean(svmlinear_mae_temp))
  svmrbf_mae = c(svmrbf_mae, mean(svmrbf_mae_temp))
  svmlinear_acc = c(svmlinear_acc, mean(svmlinear_temp))
  svmrbf_acc = c(svmrbf_acc, mean(svmrbf_temp))
}

write.csv(cbind(x, lda_acc, pda4_acc, pda6_acc, pda8_acc, svmlinear_acc, svmrbf_acc,
                   lda_mse, pda4_mse, pda6_mse, pda8_mse, svmlinear_mse, svmrbf_mse,
                   lda_mae, pda4_mae, pda6_mae, pda8_mae, svmlinear_mae, svmrbf_mae), "3_num_classes.csv")
```
## Simulation (extra): non-Gaussian data

```{r}
n = 40
p = 500
num_unimp = 490
num_imp = 10

lda_temp  = c()
pda4_temp = c()
pda6_temp = c()
pda8_temp = c()
svm_temp = c()

for (j in 1:nrun){ #compute accuracy score for 20 times
    df = scale(as.data.frame(cbind(
      rbind(matrix(rnorm(n/2*num_imp, 0, 1/4)+rnorm(n/2*num_imp, 4, 1/4),,num_imp),
            matrix(rnorm(n/2*num_imp, -4, 1/4)+rnorm(n/2*num_imp, 0, 1/4),,num_imp)),
      matrix(rnorm(n*num_unimp, -2, 1/4)+rnorm(n*num_unimp, 2, 1/4),,num_unimp))))
    class = c(rep(1, n/2), rep(2, n/2))

    ind = sample(1:40, 10, replace=FALSE)
    test = df[ind,]
    train = df[-ind,]
    cls_test = class[ind]
    cls_train = class[-ind]
  
    PP.opt = PP.optimize.anneal("LDA",2,train,cls_train)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    lda_temp = c(lda_temp, misclass(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.4)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda4_temp = c(pda4_temp, misclass(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.6)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda6_temp = c(pda6_temp, misclass(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    PP.opt = PP.optimize.anneal("PDA",2,train,cls_train,lambda = 0.8)
    proj.data.test = as.matrix(test)%*%PP.opt$proj.best
    proj.data.train = as.matrix(train)%*%PP.opt$proj.best
    pda8_temp = c(pda8_temp, misclass(PP.opt$proj.best, t(train), t(test), cls_train, cls_test))
    
    
  classifier = svm(formula = cls_train ~ .,
                 data = train,
                 type = 'C-classification',
                 kernel = 'linear')
    y_pred = predict(classifier, newdata = test)
    svm_temp = c(svm_temp, sum(y_pred == cls_test)/nrow(test))
}
```

## PLOT

```{r}
library(latex2exp)
library(tidyverse)
library(ggplot2)
library(patchwork)
cbPalette <- c("#999999", "#0072B2", "#56B4E9", "#009E73", "#D55E00","#E69F00","#F0E442","#CC79A7")
labels = c("LDA", TeX(r'(PDA $\lambda=0.4$)'),
                  TeX(r'(PDA $\lambda=0.6$)'),
                  TeX(r'(PDA $\lambda=0.8$)'),
                  "SVM(linear)", 
                  "SVM(RBF)")
```

```{r}
df = read.csv("../data/1_perc_imp_var.csv")
df %>% pivot_longer(cols=c('lda_acc', 'pda4_acc', 
                           'pda6_acc', 'pda8_acc', 
                           'svm_linear', 'svm_rbf'),
                    names_to='types',
                    values_to='value') %>%
  ggplot(aes(x = x, y = value, group = types)) +
  geom_line(aes(linetype = types, color = types))+
  geom_point(aes(color = types, shape = types))+
  scale_color_manual(values=cbPalette,
                     labels=labels)+
  scale_linetype_discrete(labels=labels)+
  scale_shape_discrete(labels=labels)+
  theme_bw()+
  xlab("Percentage of important variables")+
  ylab("Accuracy")+
  theme(legend.title=element_blank())
```  

```{r}
df = read.csv("../data/2_ratio_dim_obs.csv")
df %>% pivot_longer(cols=c('lda_acc', 'pda4_acc', 
                           'pda6_acc', 'pda8_acc', 
                           'svm_linear', 'svm_rbf'),
                    names_to='types',
                    values_to='value') %>%
  ggplot(aes(x = x/500, y = value, group = types)) +
  geom_line(aes(linetype = types, color = types))+
  geom_point(aes(color = types, shape = types))+
  scale_color_manual(values=cbPalette,
                     labels=labels)+
  scale_linetype_discrete(labels=labels)+
  scale_shape_discrete(labels=labels)+
  theme_bw()+
  xlab("Ratio between number of observations and dimensions")+
  ylab("Accuracy")+
  theme(legend.title=element_blank()) 
```

```{r}
df = read.csv("../data/3_num_classes.csv")
df %>% pivot_longer(cols=c('lda_acc', 'pda4_acc', 
                           'pda6_acc', 'pda8_acc', 
                           'svm_linear', 'svm_rbf'),
                    names_to='types',
                    values_to='value') %>%
  ggplot(aes(x = x, y = value, group = types)) +
  geom_line(aes(linetype = types, color = types))+
  geom_point(aes(color = types, shape = types))+
  scale_color_manual(values=cbPalette,
                     labels=labels)+
  scale_linetype_discrete(labels=labels)+
  scale_shape_discrete(labels=labels)+
  theme_bw()+
  xlab("Number of classes")+
  ylab("Accuracy")+
  theme(legend.title=element_blank()) 
```

```{r}
par(mfrow = c(1, 2)) 
df = read.csv("../data/4_outliers_imp.csv")
plot1 = df %>% pivot_longer(cols=c('lda_acc', 'pda4_acc', 
                           'pda6_acc', 'pda8_acc', 
                           'svm_linear', 'svm_rbf'),
                    names_to='types',
                    values_to='value') %>%
  ggplot(aes(x = x, y = value, group = types)) +
  geom_line(aes(linetype = types, color = types))+
  geom_point(aes(color = types, shape = types))+
  scale_color_manual(values=cbPalette,
                     labels=labels)+
  scale_linetype_discrete(labels=labels)+
  scale_shape_discrete(labels=labels)+
  theme_bw()+
  xlab("Number of outliers added to important variables")+
  ylab("Accuracy")+
  theme(legend.title=element_blank()) 

df = read.csv("../data/4_outliers.csv")
plot2 = df %>% pivot_longer(cols=c('lda_acc', 'pda4_acc', 
                           'pda6_acc', 'pda8_acc', 
                           'svm_linear', 'svm_rbf'),
                    names_to='types',
                    values_to='value') %>%
  ggplot(aes(x = x, y = value, group = types)) +
  geom_line(aes(linetype = types, color = types))+
  geom_point(aes(color = types, shape = types))+
  scale_color_manual(values=cbPalette,
                     labels=labels)+
  scale_linetype_discrete(labels=labels)+
  scale_shape_discrete(labels=labels)+
  theme_bw()+
  xlab("Number of outliers added to unimportant variables")+
  ylab("Accuracy")+
  theme(legend.title=element_blank()) 
plot1/plot2+ plot_layout(guides = 'collect')
```


