Je veux mieux comprendre les packages R Lars
et Glmnet
, qui sont utilisés pour résoudre le problème Lasso:
(pour Variables et échantillons, voir www.stanford.edu/~hastie/Papers/glmnet.pdf à la page 3)
Par conséquent, je les ai appliqués tous les deux sur le même jeu de données de jouets. Malheureusement, les deux méthodes ne donnent pas les mêmes solutions pour la même entrée de données. Quelqu'un at-il une idée d'où vient la différence?
J'ai obtenu les résultats comme suit: Après avoir généré des données (8 échantillons, 12 fonctionnalités, conception Toeplitz, tout centré), j'ai calculé l'ensemble du chemin Lasso en utilisant Lars. Ensuite, j'ai exécuté Glmnet en utilisant la séquence de lambdas calculée par Lars (multipliée par 0,5) et espérais obtenir la même solution, mais je ne l'ai pas fait.
On peut voir que les solutions sont similaires. Mais comment expliquer les différences? Veuillez trouver mon code ci-dessous. Il y a une question connexe ici: GLMNET ou LARS pour calculer les solutions LASSO? , mais il ne contient pas la réponse à ma question.
Installer:
# 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 avec lambda = (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
Évidemment, si les méthodes utilisent des modèles différents, vous obtiendrez des réponses différentes. La soustraction des termes d'interception ne conduit pas au modèle sans l'interception car les meilleurs coefficients d'ajustement changeront et vous ne les modifiez pas de la façon dont vous vous en approchez. Vous devez adapter le même modèle avec les deux méthodes si vous souhaitez obtenir les mêmes ou presque les mêmes réponses.
la source
Les résultats doivent être les mêmes. Le package lars utilise par défaut type = "lar", changez cette valeur en type = "lasso". Il suffit de baisser le paramètre 'thresh = 1e-16' pour glmnet car la descente de coordonnées est basée sur la convergence.
la source