Воспроизведите пример моделирования лассо 1

Question

Я хотел бы воспроизвести результаты примера 1 на странице 280 в оригинальной статье лассо.

• Модель y = X*beta + sigma*epsilon где epsilon является N(0,1)
• Смоделируйте 50 наборов данных, состоящих из 20/20/200 наблюдений для обучающих / проверочных / тестовых наборов.
• True beta = (3, 1.5, 0, 0, 2, 0, 0, 0)
• sigma = 3

• Парная корреляция между x_i а также x_j должны быть corr(i,j) = 0.5^|i-j|

  Я использовал обучение, валидацию, тестовый подход, чтобы найти оценки test MSE, Я пытался вычислить несколько test MSE оценки вручную, чтобы проверить, нахожусь ли я на правильном пути перед повторениями моделирования. Но похоже test MSE оценки, которые я нахожу (между [9, 15]), намного больше, чем оценки, приведенные в оригинальной статье (со средним значением 2,43). Я прилагаю код, который я использовал для генерации test MSE"S.

Любое предложение, пожалуйста?

    library(MASS)
    library(glmnet)

    simfun <- function(trainn = 20, validationn = 20, testn = 200, corr =0.5, sigma = 3, beta) { 


      n <- trainn + testn + validationn
      p <- length(beta)
      Covmatrix <- outer(1:p, 1:p, function(x,y){corr^abs(x-y)})
      X <- mvrnorm(n, rep(0,p), Covmatrix) # MASS
      X <- X[sample(n),]
      y <- X%*%beta + rnorm(n,mean = 0,sd=sigma)
      trainx <- X[1:trainn,]
      validationx <- X[(trainn+1):(trainn+validationn),]
      testx <- X[(trainn+validationn+1):n,]
      trainy <- y[1:trainn,]
      validationy <- y[(trainn+1):(trainn+validationn),]
      testy <- y[(trainn+validationn+1):n,]
      list(trainx = trainx, validationx = validationx, testx = testx, 
           trainy = trainy, validationy = validationy, testy = testy)
    }

    beta <- c(3,1.5,0,0,2,0,0,0)
    data <- simfun(20,20,200,corr=0.5,sigma=3,beta)
    trainx <- data$trainx
    trainy <- data$trainy
    validationx <- data$validationx
    validationy <- data$validationy
    testx <- data$testx
    testy <- data$testy


    # training: find betas for all the lambdas
    betas <- coef(glmnet(trainx,trainy,alpha=1))

    # validation: compute validation test error for each lambda and find the minimum
    J.val <- colMeans((validationy-cbind(1,validationx)%*%betas)^2)
    beta.opt <- betas[, which.min(J.val)]

    # test
    test.mse <- mean((testy-cbind(1,testx)%*%beta.opt)^2)
    test.mse

Answer 1

Это симуляционное исследование, поэтому я думаю, что вам не нужно использовать подход обучения-проверки. Это просто вызывает вариации из-за своей случайности. Вы можете реализовать ожидаемую ошибку теста, используя ее определение.

1. Создайте несколько наборов тренировочных данных после вашей конструкции
2. Создать независимый набор тестов
3. Подгонка каждой модели на основе каждого тренировочного набора
4. Вычислить погрешность по сравнению с тестовым набором

5. Возьми среднее

   set.seed(1)
   simpfun <- function(n_train = 20, n_test = 10, simul = 50, corr = .5, sigma = 3, beta = c(3, 1.5, 0, 0, 2, 0, 0, 0), lam_grid = 10^seq(-3, 5)) {
     require(foreach)
     require(tidyverse)
     # true model
     p <- length(beta)
     Covmatrix <- outer(
       1:p, 1:p,
       function(x, y) {
         corr^abs(x - y)
       }
     )
     X <- foreach(i = 1:simul, .combine = rbind) %do% {
       MASS::mvrnorm(n_train, rep(0, p), Covmatrix)
     }
     eps <- rnorm(n_train, mean = 0, sd = sigma)
     y <- X %*% beta + eps # generate true model
     # generate test set
     test <- MASS::mvrnorm(n_test, rep(0, p), Covmatrix)
     te_y <- test %*% beta + rnorm(n_test, mean = 0, sd = sigma) # test y
     simul_id <- gl(simul, k = n_train, labels = 1:n_train)
     # expected test error
     train <-
       y %>%
       as_tibble() %>%
       mutate(m_id = simul_id) %>%
       group_by(m_id) %>% # for each simulation
       do(yhat = predict(glmnet::cv.glmnet(X, y, alpha = 1, lambda = lam_grid), newx = test, s = "lambda.min")) # choose the smallest lambda
     MSE <- # (y0 - fhat0)^2
       sapply(train$yhat, function(x) {
         mean((x - te_y)^2)
       })
     mean(MSE) # 1/simul * MSE
   }
   simpfun()
   

---

Редактировать: для настройки параметра,

    find_lambda <- function(.data, x, y, lambda, x_val, y_val) {
      .data %>%
        do(
          tuning = predict(glmnet::glmnet(x, y, alpha = 1, lambda = lambda), newx = x_val)
        ) %>%
        do( # tuning parameter: validation set
          mse = apply(.$tuning, 2, function(yhat, y) {
            mean((y - yhat)^2)
          }, y = y_val)
        ) %>%
        mutate(mse_min = min(mse)) %>%
        mutate(lam_choose = lambda[mse == mse_min]) # minimize mse
    }

Используя эту функцию, можно добавить шаг проверки

    simpfun <- function(n_train = 20, n_val = 20, n_test = 10, simul = 50, corr = .5, sigma = 3, beta = c(3, 1.5, 0, 0, 2, 0, 0, 0), lam_grid = 10^seq(10, -1, length.out = 100)) {
    require(foreach)
    require(tidyverse)
    # true model
    p <- length(beta)
    Covmatrix <- outer(
      1:p, 1:p,
      function(x, y) {
        corr^abs(x - y)
      }
    )
    X <- foreach(i = 1:simul, .combine = rbind) %do% {
      MASS::mvrnorm(n_train, rep(0, p), Covmatrix)
    }
    eps <- rnorm(n_train, mean = 0, sd = sigma)
    y <- X %*% beta + eps # generate true model
    # generate validation set
    val <- MASS::mvrnorm(n_val, rep(0, p), Covmatrix)
    val_y <- val %*% beta + rnorm(n_val, mean = 0, sd = sigma) # validation y
    # generate test set
    test <- MASS::mvrnorm(n_test, rep(0, p), Covmatrix)
    te_y <- test %*% beta + rnorm(n_test, mean = 0, sd = sigma) # test y
    simul_id <- gl(simul, k = n_train, labels = 1:n_train)

    Y <-
      y %>%
      as_tibble() %>%
      mutate(m_id = simul_id) %>%
      group_by(m_id) %>% # for each simulation: repeat
      rename(y = V1)

    # Tuning parameter
    Tuning <-
      Y %>%
      find_lambda(x = X, y = y, lambda = lam_grid, x_val = val, y_val = val_y)

    # expected test error
    test_mse <-
      Tuning %>%
      mutate(
        test_err = mean(
          (predict(glmnet::glmnet(X, y, alpha = 1, lambda = lam_choose), newx = test) - te_y)^2
        )
      ) %>%
      select(test_err) %>%
      pull()
    mean(test_mse)
  }
  simpfun()