0

我在模拟数据集上拟合了一个模型来比较 glmnet 和 CVXR 结果。

如果我没有代码错误，结果就会大不相同。

明确地 glmnet 产生非常接近真实参数的结果。

为什么会这样？

library(CVXR)
library(glmnet)

set.seed(571)
n = 500
p = 9
x = matrix(rnorm(n*p), ncol=p)
b = c(0.5, 0, 25, -25, 125, -125, rep(0, 3))
y = x %*% b + rnorm(n, sd=.05)

n = nrow(x); p = ncol(x)

lam = 0.4
al = 0.3

# glmnet

glmnet_res = coef(glmnet(x,y,alpha=al,standardize=F,intercept=F),s=lam)[-1]

# CVXR

elastic_reg = function(beta, lambda = 0, alpha = 0) {
  ridge = 0.5*(1 - alpha) * sum(beta^2)
  lasso = alpha * p_norm(beta, 1)
  lambda * (lasso + ridge)
}

beta = Variable(p)  
loss = sum((y - x %*% beta)^2)/(2*n)

## Elastic-net regression
obj = loss + elastic_reg(beta, lam, al)
prob = Problem(Minimize(obj))
result = solve(prob)
beta_vals = result$getValue(beta)

cvxr_res = round(beta_vals,7)

cbind(glmnet_res,cvxr_res)

结果

      glmnet_res    cvxr_res         
 [1,]    0.00000   0.2417734
 [2,]    0.00000   0.0000475
 [3,]   23.39102  19.0372445
 [4,]  -23.26282 -18.6020795
 [5,]  121.59156  96.7286536
 [6,] -121.17658 -95.0466518
 [7,]    0.00000  -1.8589296
 [8,]    0.00000   0.2651426
 [9,]    0.00000   1.0167725
4

1 回答 1

0

对于连续结果，glmnet 通过其标准差缩放结果 (y)。将 glmnet 中的解决方案与其他软件进行比较的最简单方法是显式缩放 y。此外，您需要lam按标准差缩放您在 CVXR 中使用的相应惩罚值 ( )，因为您提供的惩罚值coef()也会自动按 y 的标准差缩放。在拟合之后，来自 CVXR 的估计值可以是非标准化的。我还对您的代码做了另外两个小改动：

  • 更改了glmnet和 CVXR 的收敛阈值以提高准确性
  • 增加了惩罚值lam（

修改后的代码

library(CVXR)
library(glmnet)

# simulate data
set.seed(571)
n <- 500
p <- 9
x <- matrix(rnorm(n*p), ncol=p)
b <- c(0.5, 0, 25, -25, 125, -125, rep(0, 3))
y <- x %*% b + rnorm(n, sd = .5)
sd_y <- drop(sqrt(var(y) * (n - 1) / n))
y_stnd <- y / sd_y

# fix penalty value and EN parameter
lam <- 20
al <- 0.3

# fit EN in glmnet
fit_glmnet <- glmnet(x = x,
                     y = y,
                     alpha = al,
                     standardize = FALSE, 
                     intercept = FALSE,
                     thresh = 1e-20)

betas_glmnet <- as.vector(coef(fit_glmnet, 
                               s = lam, 
                               exact = TRUE, 
                               x = x, 
                               y = y)[-1])

# fit EN in CVXR (using standardized y and rescaled penalty, lambda / sd_y)
beta <- Variable(p)  
obj <- Minimize(sum((y_stnd - x %*% beta)^2) / (2 * n) + 
                (lam / sd_y) * ((1 - al) * sum_squares(beta) / 2 + al * p_norm(beta, 1)))
prob <- Problem(obj)
result <- solve(prob, solver = "ECOS", verbose = TRUE, ABSTOL = 1e-12, RELTOL = 1e-10)
betas_cvxr <- drop(result$getValue(beta))

# Compare results (unstandardize estimates for CVXR)
round(cbind(betas_glmnet, sd_y * betas_cvxr), 6)

结果

 [1,]      0.00000    0.00000
 [2,]      0.00000    0.00000
 [3,]     17.84706   17.84706
 [4,]    -17.28221  -17.28221
 [5,]    109.82539  109.82539
 [6,]   -108.07262 -108.07262
 [7,]      0.00000    0.00000
 [8,]      0.00000    0.00000
 [9,]      0.00000    0.00000
于 2018-11-16T21:05:19.617 回答