Chcę lepiej zrozumieć pakiety R Lars
i Glmnet
używane do rozwiązania problemu Lasso:
(dla zmiennych i próbek , patrz www.stanford.edu/~hastie/Papers/glmnet.pdf na stronie 3)
Dlatego zastosowałem je oba w tym samym zestawie danych zabawek. Niestety obie metody nie dają takich samych rozwiązań dla tych samych danych wejściowych. Czy ktoś ma pomysł, skąd ta różnica?
Otrzymałem wyniki w następujący sposób: Po wygenerowaniu niektórych danych (8 próbek, 12 cech, projekt Toeplitz, wszystko wyśrodkowane), obliczyłem całą ścieżkę Lasso za pomocą Larsa. Następnie uruchomiłem Glmnet, używając sekwencji lambd obliczonej przez Larsa (pomnożonej przez 0,5) i liczyłem na uzyskanie tego samego rozwiązania, ale nie zrobiłem tego.
Widać, że rozwiązania są podobne. Ale jak mogę wyjaśnić różnice? Mój kod znajduje się poniżej. Tutaj jest powiązane pytanie: GLMNET lub LARS do obliczania rozwiązań LASSO? , ale nie zawiera odpowiedzi na moje pytanie.
Ustawiać:
# 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 z 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
Oczywiście, jeśli metody wykorzystują różne modele, otrzymasz różne odpowiedzi. Odjęcie warunków przechwytywania nie prowadzi do modelu bez przechwytywania, ponieważ zmienią się współczynniki najlepszego dopasowania i nie zmienisz ich w miarę zbliżania się do niego. Musisz dopasować ten sam model do obu metod, jeśli chcesz uzyskać te same lub prawie takie same odpowiedzi.
źródło
Wyniki muszą być takie same. pakiet lars używa domyślnie type = "lar", zmień tę wartość na type = "lasso". Po prostu obniż parametr „thresh = 1e-16” dla glmnet, ponieważ opadanie współrzędnych opiera się na zbieżności.
źródło