0
$\begingroup$

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 yi be equally divided between classes 0 and 1, and let xi∈R10 be normally distributed.

Given yi=0, xi∼N10(0,I10). Given yi=1, xi∼N10(μ,I10) with μ=(1,1,1,1,1,0,0,0,0,0).

The N10 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.

$\endgroup$
0
$\begingroup$

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 .1634

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.