Я хочу лучше понять пакеты R Lars
и Glmnet
, которые используются для решения проблемы Лассо:
(для переменных и выборок, см. www.stanford.edu/~hastie/Papers/glmnet.pdf на стр. 3)
Поэтому я применил их обоих к одному и тому же набору игрушечных данных. К сожалению, оба метода не дают одинаковых решений для одного и того же ввода данных. У кого-нибудь есть идея, откуда берется эта разница?
Я получил результаты следующим образом: После генерации некоторых данных (8 образцов, 12 объектов, дизайн Тёплица, все по центру) я вычислил весь путь Лассо, используя Lars. Затем я запустил Glmnet, используя последовательность лямбда-выражений, вычисленную Ларсом (умноженную на 0,5), и надеялся получить то же решение, но я этого не сделал.
Видно, что решения похожи. Но как я могу объяснить различия? Пожалуйста, найдите мой код ниже. Здесь есть связанный вопрос: GLMNET или LARS для вычисления решений LASSO? , но он не содержит ответа на мой вопрос.
Настроить:
# Load packages.
library(lars)
library(glmnet)
library(MASS)
# Set parameters.
nb.features <- 12
nb.samples <- 8
nb.relevant.indices <- 3
snr <- 1
nb.lambdas <- 10
# Create data, not really important.
sigma <- matrix(0, nb.features, nb.features)
for (i in (1:nb.features)) {
for (j in (1:nb.features)) {
sigma[i, j] <- 0.99 ^ (abs(i - j))
}
}
x <- mvrnorm(n=nb.samples, rep(0, nb.features), sigma, tol=1e-6, empirical=FALSE)
relevant.indices <- sample(1:nb.features, nb.relevant.indices, replace=FALSE)
x <- scale(x)
beta <- rep(0, times=nb.features)
beta[relevant.indices] <- runif(nb.relevant.indices, 0, 1)
epsilon <- matrix(rnorm(nb.samples),nb.samples, 1)
simulated.snr <-(norm(x %*% beta, type="F")) / (norm(epsilon, type="F"))
epsilon <- epsilon * (simulated.snr / snr)
y <- x %*% beta + epsilon
y <- scale(y)
Lars:
la <- lars(x, y, intercept=TRUE, max.steps=1000, use.Gram=FALSE)
co.lars <- as.matrix(coef(la, mode="lambda"))
print(round(co.lars, 4))
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] 0.0000 0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
# [2,] 0.0000 0 0 0.0000 0.0000 0.1735 0.0000 0.0000 0.0000 0.0000
# [3,] 0.0000 0 0 0.2503 0.0000 0.4238 0.0000 0.0000 0.0000 0.0000
# [4,] 0.0000 0 0 0.1383 0.0000 0.7578 0.0000 0.0000 0.0000 0.0000
# [5,] -0.1175 0 0 0.2532 0.0000 0.8506 0.0000 0.0000 0.0000 0.0000
# [6,] -0.3502 0 0 0.2676 0.3068 0.9935 0.0000 0.0000 0.0000 0.0000
# [7,] -0.4579 0 0 0.6270 0.0000 0.9436 0.0000 0.0000 0.0000 0.0000
# [8,] -0.7848 0 0 0.9970 0.0000 0.9856 0.0000 0.0000 0.0000 0.0000
# [9,] -0.3175 0 0 0.0000 0.0000 3.4488 0.0000 0.0000 -2.1714 0.0000
# [10,] -0.4842 0 0 0.0000 0.0000 4.7731 0.0000 0.0000 -3.4102 0.0000
# [11,] -0.4685 0 0 0.0000 0.0000 4.7958 0.0000 0.1191 -3.6243 0.0000
# [12,] -0.4364 0 0 0.0000 0.0000 5.0424 0.0000 0.3007 -4.0694 -0.4903
# [13,] -0.4373 0 0 0.0000 0.0000 5.0535 0.0000 0.3213 -4.1012 -0.4996
# [14,] -0.4525 0 0 0.0000 0.0000 5.6876 -1.5467 1.5095 -4.7207 0.0000
# [15,] -0.4593 0 0 0.0000 0.0000 5.7355 -1.6242 1.5684 -4.7440 0.0000
# [16,] -0.4490 0 0 0.0000 0.0000 5.8601 -1.8485 1.7767 -4.9291 0.0000
# [,11] [,12]
# [1,] 0.0000 0.0000
# [2,] 0.0000 0.0000
# [3,] 0.0000 0.0000
# [4,] -0.2279 0.0000
# [5,] -0.3266 0.0000
# [6,] -0.5791 0.0000
# [7,] -0.6724 0.2001
# [8,] -1.0207 0.4462
# [9,] -0.4912 0.1635
# [10,] -0.5562 0.2958
# [11,] -0.5267 0.3274
# [12,] 0.0000 0.2858
# [13,] 0.0000 0.2964
# [14,] 0.0000 0.1570
# [15,] 0.0000 0.1571
glmnet с лямбда = (lambda_lars / 2):
glm2 <- glmnet(x, y, family="gaussian", lambda=(0.5 * la$lambda), thresh=1e-16)
co.glm2 <- as.matrix(t(coef(glm2, mode="lambda")))
print(round(co.glm2, 4))
# (Intercept) V1 V2 V3 V4 V5 V6 V7 V8 V9
# s0 0 0.0000 0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
# s1 0 0.0000 0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
# s2 0 0.0000 0 0 0.2385 0.0000 0.4120 0.0000 0.0000 0.0000
# s3 0 0.0000 0 0 0.2441 0.0000 0.4176 0.0000 0.0000 0.0000
# s4 0 0.0000 0 0 0.2466 0.0000 0.4200 0.0000 0.0000 0.0000
# s5 0 0.0000 0 0 0.2275 0.0000 0.4919 0.0000 0.0000 0.0000
# s6 0 0.0000 0 0 0.1868 0.0000 0.6132 0.0000 0.0000 0.0000
# s7 0 -0.2651 0 0 0.2623 0.1946 0.9413 0.0000 0.0000 0.0000
# s8 0 -0.6609 0 0 0.7328 0.0000 1.6384 0.0000 0.0000 -0.5755
# s9 0 -0.4633 0 0 0.0000 0.0000 4.6069 0.0000 0.0000 -3.2547
# s10 0 -0.4819 0 0 0.0000 0.0000 4.7546 0.0000 0.0000 -3.3929
# s11 0 -0.4767 0 0 0.0000 0.0000 4.7839 0.0000 0.0567 -3.5122
# s12 0 -0.4715 0 0 0.0000 0.0000 4.7915 0.0000 0.0965 -3.5836
# s13 0 -0.4510 0 0 0.0000 0.0000 5.6237 -1.3909 1.3898 -4.6583
# s14 0 -0.4552 0 0 0.0000 0.0000 5.7064 -1.5771 1.5326 -4.7298
# V10 V11 V12
# s0 0.0000 0.0000 0.0000
# s1 0.0000 0.0000 0.0000
# s2 0.0000 0.0000 0.0000
# s3 0.0000 0.0000 0.0000
# s4 0.0000 0.0000 0.0000
# s5 0.0000 -0.0464 0.0000
# s6 0.0000 -0.1293 0.0000
# s7 0.0000 -0.4868 0.0000
# s8 0.0000 -0.8803 0.3712
# s9 0.0000 -0.5481 0.2792
# s10 0.0000 -0.5553 0.2939
# s11 0.0000 -0.5422 0.3108
# s12 0.0000 -0.5323 0.3214
# s13 -0.0503 0.0000 0.1711
# s14 0.0000 0.0000 0.1571
Очевидно, что если методы используют разные модели, вы получите разные ответы. Вычитание членов перехвата не приводит к модели без перехвата, потому что изменятся наилучшие подходящие коэффициенты, и вы не измените их так, как подходите к ним. Вам нужно использовать одну и ту же модель для обоих методов, если вы хотите получить одинаковые или почти одинаковые ответы.
источник
Результаты должны быть одинаковыми. Пакет lars по умолчанию использует type = "lar", измените это значение на type = "lasso". Просто уменьшите параметр 'thresh = 1e-16' для glmnet, так как спуск по координатам основан на сходимости.
источник