R语言实现LASSO回归——自己编写LASSO回归算法

这篇文章中我们可以编写自己的代码来计算套索(lasso)回归

由Kaizong Ye,Coin Ge撰写

我们必须定义阈值函数

×

(1)尽量减少选取变量的数量

可以人工检查每一项变量,并以此来确定哪些变量更为重要,保留那些更为重要的特征变量。这种做法非常有效,但是其缺点是当舍弃一部分特征变量时,也舍弃了问题中的一些信息。例如,所有的特征变量对于预测房价都是有用的,实际上并不想舍弃一些信息或者说舍弃这些特征变量。

(2)正则化

在正则化中将保留所有的特征变量,但是会减小特征变量的数量级。 当有很多特征变量时,其中每一个变量都能对预测产生一点影响。假设在房价预测的例子中,可以有很多特征变量,其中每一个变量都是有用的,因此不希望把任何一个变量删掉,这就导致了正则化概念的发生。

如果应用正则化去解决实际问题?

在图1中,可以看出最右边的是过拟合的,但是不想抛弃\theta_{3}\theta_{4}变量自带的信息,因此加上惩罚项,使得\theta_{3}\theta_{4}足够小。

优化目标为:

                                                     \min _{\theta} \frac{1}{2 m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right)^{2}

在优化目标的基础上,添加一些项(惩罚项),如下:

                                      \min _{\theta} \frac{1}{2 m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right)^{2}+1000 \theta_{3}^{2}+1000 \theta_{4}^{2}

1000 只是随便写的一个表示很大的数字。现在,如果想要最小化这个函数,那么为了最小化这个新的损失函数,就要让 \theta_{3}\theta_{4}尽可能小。因为,如果在原有代价函数的基础上加上 1000 \theta_{3}^{2}+1000 \theta_{4}^{2}这一项 ,那么这个新的代价函数将变得很大,所以,当最小化这个新的代价函数时, \theta_{3}\theta_{4}尽可能小或者接近于0,就像我们忽略了这两个值一样。如果做到这一点(\theta_{3}\theta_{4}接近于0),那么我们将得到一个近似的二次函数,并且h_{\theta}中并没有丢失\theta_{3}\theta_{4}的信息(虽然\theta_{3}\theta_{4}接近于0,但是还是存在的)

但是在真实问题中,一个预测问题可能有上百种特征,并不知道哪一些是高阶多项式的项,就像不知道\theta_{3}\theta_{4}是高阶的项。所以,如果假设有一百个特征,并不知道如何选择关联度更好的参数,如何缩小参数的数目等等。 因此在正则化里,要做的事情,就是把减小损失函数中所有的参数值,(因为并不知道是哪一个或哪几个是要去缩小的)。

因此,我们需要修改代价函数,变成如下形式:

 

                                            \frac{1}{2 m}\left[\sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right)^{2}+\lambda \sum_{j=1}^{n} \theta_{j}^{2}\right]   

其中\lambda \sum_{j=1}^{n} \theta_{j}^{2}就是正则化项; λ称为正则化参数,作用是控平衡拟合训练的目标和保持参数值较小。

Lasso回归模型

Lasso回归是在损失函数后,加L1正则化,如下所示:

                                              \frac{1}{2 m}\left[\sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right)^{2}+ \lambda \sum_{j=1}^{k}\left|w_{j}\right|\right]

m为样本个数,k为参数个数,其中\lambda \sum_{j=1}^{k}\left|w_{j}|\right为L1正则化。

此外:还有L2正则化:\lambda \sum_{j=1}^{n} w_{j}^{2},加了这种正则化的损失函数模型为:脊(岭)回归(Ridge Regression)。



R函数是

thresh = function(x,a){
sign(x) * pmax(abs(x)-a,0)
}

要解决我们的优化问题,设置


视频

Lasso回归、岭回归等正则化回归数学原理及R语言实例

探索见解

去bilibili观看

探索更多视频

这样就可以等效地写出优化问题

因此

一个得到

同样,如果有权重ω=(ωi),则按坐标更新将变为

计算此分量下降的代码是

lasso = function(X,y,beta,lambda,tol=1e-6,maxiter=1000){
 
beta0 = sum(y-X%*%beta /(length(y))
beta0list[1] = beta0
for (j in 1:maxiter){
for (k in 1:length beta)){
r = y - X[,-k]%*%beta[-k] - beta0*rep(1,length(y )
beta[k] = (1/sum(omega*X[,k]^2) *
threshog(t(omega*r)%*%X[,k ,length(y *lambda)
}
beta0 = sum(y-X%*%beta)/(length(y))
 
 
obj[j] = (1/2)*(1/length(y))*norm(omega*(y - X%*%beta - 
beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta))
if (norm(rbind(beta0list[j],betalist[[j]]) - 
rbind(beta0,beta),'F') ) { break } 

 例如,考虑以下(简单)数据集,其中包含三个协变量


chicago = read.table("data.txt",header=TRUE,sep=";")

我们可以“标准化”

 
for(j in 1:3) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j])
 
y = (y-mean(y))/sd(y)

 要初始化算法,使用OLS估算

lm(y~0+.,)$coef

例如

 
lasso(X,y,beta_init,lambda=.001)
$obj
[1] 0.001014426 0.001008009 0.001009558 0.001011094 0.001011119 0.001011119
 
$beta
[,1]
X_1 0.0000000
X_2 0.3836087
X_3 -0.5026137
 
$intercept
[1] 2.060999e-16

 我们可以通过循环获得标准的lasso图


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498