r - R 中的这个套索回归实现有什么问题?
问题描述
我正在尝试从头开始编写一个函数,该函数使用具有协调下降(Gauss-Seidel)和软阈值的 LASSO 估计回归系数。
我的代码如下:
library(MASS)
set.seed(1)
n = 200
p = 200
V = matrix(0.2, p, p)
diag(V) = 1
X = as.matrix(mvrnorm(n, mu = rep(0, p), Sigma = V))
y = X[, 1] + 0.5*X[, 2] + 0.25*X[, 3] + rnorm(n)
X = scale(X)
y = scale(y)
soft_th <- function(b, lambda){
if (b > lambda){
return (b - lambda)
}
else if (b < -lambda){
return(b + lambda)
}
else {
return (0)
}
}
myLasso <- function(X,y, lambda=0.3,tol=1e-5,maxitr=100){
beta_old <- rep(0,p)
beta_new <- rep(0,p)
for(i in 1:maxitr){
beta_old <- beta_new
for (j in (1:ncol(X)))
{
X_j <- X[,j]
y_pred <- t(X)%*%beta_old
rho <- t(X_j)%*%(y - y_pred + beta_old[j]*X_j)
beta_new[j] <- soft_th(rho,0.7)
}
l1 <- sum(abs(beta_old-beta_new))
print(l1)
r <- y - t(X)%*%beta_old
if (l1<tol){
print('Convergence reached')
break
}
}
}
myLasso(X,y)
我遇到的问题是 beta_old 和 beta_new 之间的 L1 范数在每次迭代之间增加(很多!)。我正在关注这篇文章中所说的内容:
我认为我没有正确实施下降更新规则。
任何帮助,将不胜感激。提前致谢。
解决方案
我正在做更多的研究,似乎我没有对 X 矩阵进行归一化。定义X后加上X <- X/norm(X,type='2'),问题就解决了。
现在我遇到的一个新问题是,这个新函数不会复制 glmnet 实现 LASSO 回归的结果。会是什么呢?glmnet 的 RMSE 为 0.6,我的实现为 0.997。如果有人可以指导我如何改进我的功能,我会很高兴。
提前致谢。
推荐阅读
- python - 使用张量流进行 CNN 预测
- javascript - React-DND 不会触发下降并悬停在目标上
- azure - 功能级别授权 Azure 中的授权密钥我们可以通过 CICD 管理这些密钥吗
- c++ - C++:对 ABC 的子类强制执行单例
- r - 如何安全地将 api 密钥添加到 dockerized Shiny App?
- python - Python Django 错误 ModuleNotFoundError:没有名为“corsheaders”的模块
- python - 为什么 np.ones 创建黑色图像?
- python - 在打印的字符串中打印希腊字母
- c - Eclipse CDT:如何从源文件生成头文件?
- ruby - Bundler:将 gem 从 git 安装到 gems 文件夹而不是 bundler 文件夹