SVM and Monte Carlo simulation to compute misclassification error rate

Question

I am trying to solve the following problem with R: use simulation to evaluate (by Monte Carlo) the expected misclassification error rate given a particular generating model. Let $y_i$ be equally divided between classes $0$ and $1$, and let $x_i \in \mathbb{R}^{10}$ be normally distributed.

Given $y_i=0$, $x_i ∼ N_{10}(0, I_{10})$. Given $y_i=1$, $x_i ∼ N_{10}(\mu, I_{10})$ with $\mu = (1, 1, 1, 1, 1, 0, 0, 0, 0, 0)$.

The $N_{10}$ notation just means its a ten-dimensional Gaussian distribution; you can use the mvrnorm function in the MASS package to help generate the data. Now, we would like to know the expected test error rate if we fit an SVM to a sample of 50 random training points from class 1 and 50 more from class 0. We can calculate this to high precision by 1) generating a random training sample to train on, 2) evaluating the number of mistakes we make on a large test set, and then 3) repeating (1-2) many times and averaging the error rate for each trial.

Use svm in the e1071 package with the default settings (the default kernel is a radial kernel). What is the expected test error rate of this method.

I used two different approach which give similar results (around 0.07), but unfortunately neither of them seem to work:

    set.seed(2077)
    library(MASS)
    library(e1071)
    
    mse = rep(NA, 100)
    for (i in 1:100){
      # x = matrix(rnorm(100 * 10), ncol = 10)
      x = mvrnorm(100, mu=c(1,1,1,1,1,0,0,0,0,0), Sigma=diag(10))
      y = c(rep(1,50), rep(0,50))
      x[y==1,] = x[y==1,] + 1
      df = data.frame(x, y=as.factor(y))
      
      svm.fit = svm(y~., data=df, kernel='radial')
      
      x.test = mvrnorm(100, mu=c(1,1,1,1,1,0,0,0,0,0), Sigma=diag(10))
      # x.test = matrix(rnorm(100 * 10), ncol = 10)
      y.test = sample(c(0,1), 100, replace=T)
      x.test[y.test==1,] = x.test[y.test==1,] + 1
      df.test = data.frame(x.test, y.test=as.factor(y.test))
      svm.pred = predict(svm.fit, df.test)
      tab = table(svm.pred, df.test$y)
      mse[i] = (tab[2] + tab[3])/100
    }
    mean(mse)

    library(MASS)
    library(e1071)
    # generating the distribution
    mean = c(1,1,1,1,1,0,0,0,0,0)
    x = mvrnorm(200, mean, Sigma=diag(length(mean)))
    y = c(rep(0, 100), rep(1,100))
    x[y==1,] = x[y==1,] + 1
    df = data.frame(x=x, y=as.factor(y))
    
    # finding the seeds that split train and test properly
    seed0 = rep(NA, 1000)
    for (i in 1:1000){
      set.seed(i)
      samp = sample(1:200, 100)
      df.train = df[samp,]
      if (table(df.train$y)[1]==50){
        seed0[i] = i
      }
    }
    seed0 = na.omit(seed0)
    seed0 = seed0[1:100]
    
    # computing the expected test MSE
    mse.vector = rep(NA, 100)
    for (i in 1:100){
      set.seed(seed0[i])
      train = sample(1:200, 100)
      df.train = df[train,]
      df.test = df[-train,]
      svmfit = svm(y~., data=df.train, kernel='radial')
      svm.pred = predict(svmfit, newdata=df.test)
      tab = table(svm.pred, df.test$y)
      mse = (tab[2] + tab[3])/100
      mse.vector[i] = mse
    }
    mean(mse.vector)

Since it's a quiz I don't have the answer, but I found another post on the same question and the person who posted it says that it should be around 0.16. I hope you might help me. Thanks for your time.

Answer 1

I think your first approach is on the right track. There are two changes:

1. The distribution means are slightly different for class 0 and class 1, therefore, we need to generate them separately.

2. Instead of re-sampling from the dataset, we can generate the training set (n = 100) and test set (n = 1000 or a large number you like) separately because we know the underlying distributions in this case.

The codes below should give you an average mse = 0.164075.

library(MASS)
library(e1071)

set.seed(123)

# custom function to quickly generate dataset 
create_dataset <- function(n = 50) {
  mu <- rep(c(1, 0), c(5, 5))
  

  x0 <-
    matrix(rnorm(n * 10), ncol = 10)
  
  x1 <- 
    mvrnorm(n = n, mu = mu, Sigma = diag(10))
  
  x <- rbind(x1, x0)
  
  y <- 
    c(rep(1, n), rep(0, n))
  
  df = data.frame(x, y=as.factor(y))
  
  return(df)
}

# create vector to capture mean standard error
mse <-  rep(NA, 100)


for (i in 1:100) {
  # create 100 sample with 50 class 1 and 50 class 0
  df_train <- create_dataset(50)

  # create a large test dataset
  df_test <- create_dataset(1000)
  
  svm.fit = svm(y~., data=df_train, kernel = 'radial')
  
  svm.pred <- predict(svm.fit, newdata = df_test)
  
  tab <- table(svm.pred, df_test$y)
  
  mse[i] <- (tab[2] +tab[3])/sum(tab)
  
  
}

mean(mse)

Answer 2

On your first approach, change the line:

x[y==1,] = x[y==1,] + 1

into

x[y==1, 1:5] = x[y==1, 1:5] + 1

to allow only the first 5 columns to have a class one on your training data.

You should get an MSE around 0.1634.