Воспроизведите пример моделирования лассо 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.

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

    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()