Lasso算法与AIC、bic

Lasso算法与AIC、BIC、Stepwise算法比较

一、变量选择

回归分析中如果有很多个变量，但不进行变量选择，会使回归系数的精度下降，模型的准确率降低，还会造成统计研究的成本较大。所以变量选择在回归分析中是一个重点问题。

在回归方程中，预测精度和可解释性是评估回归模型的两个重要指标。传统的变量选择方法有forward法，toward法，逐步回归法，全子集法。结合的变量选择的标准有AIC、BIC、CP准则等。

Robert提出的Lasso回归是一种收缩估计方法，基本思想是在一个回归系数的绝对值之和小于一个常数的条件下，使残差平方和最小化，从而较严格地使系数本该为零化为零，相应的变量被删除，实现变量选择。

二、模拟实验

生成200个独立同分布的随机变量X，服从Xi~N(0,I6)，残差εi~N 0,4 ，

令β=[1,0.8,0.6,0,0,0]。然后由Yi=Xiβ+εi得到因变量Y，现已知自变量X和因变量Y，要求做变量选择， 重复模拟实验1000次，统计1000次结果中β1，β2，β3三个系数为零的个数和的平均值，记为incorrect,它是指惩罚过重将把不该估为零的系数估为零的可能性，即错误率； β4，β5，β6三个系数为零的个数和的平均值，记为correct，它是指通过惩罚函数把该估为零的系数估为零的可能性，即正确率。

最后比较四种算法的incorrect和correct.

1、LASSO算法的编程实现：

2

LASSO的选择标准是使n ni yi?xiβi +λ j|βj|达到最小。

1

p

为满足上面这一标准，现在用得较普遍的是LASSO的lars算法，在本报告我们采用LASSO中的LQA算法，主要体现了迭代的思想。迭代的步骤是：

第一步：设参数 β估计的初值为β0 ，固定λ ，选取最小二乘估计 βLS作为 β的初值，

第二步：由 估计的初值为 和下面的迭代公式得出固定 时的 ，记为 βλ

β（1）=β（0）?[?2QLSβ 0 +n β(0)]?1[?QLSβ 0 +nUλβ(0)]

λ

其中，

?2QLSβ 0 =X′X ?QLSβ 0 =?X(Y?Xβ(0)) 11

β 0 =diag ，…，

β0 βp λ

如果βj=0，那么取 β 为某个很大的值

j

1

Uλβ(0)= β 0 β(0)

λ

在我的程序中，当|βj|<0.0001时，令 β =100000

j

1

第三步：如果某个βj很小时，直接令βj=0.

在我的程序中，经过反复试验，发现取小于0.1时，另βj=0效果理想。

重复上述第二步和第三步，直到β收敛，即|β(n+1)?β(n)|<0.000001,这时得到固定一个λ的βλ估计

第四步：使λ在0~1之间变化，求出使广义交叉验证统计量

|Y?Xβλ|2

GCV λ = 1?e(λ)/n 2最小的λ，这时的λ为最优选择的罚参数，记λ?=argminGCV(λ),从而参数β的Lasso估计βlasso=βλ?

在我的程序里，经多次试验得出，λ?一般在0.3左右，为进一步细化λ，则使λ在0.2~0.5之间变化，取其步长为0.05.然后用第四步的思想求出βlasso 2、AIC和BIC准则：

AIC准则：选择子集A，使AIC A =ln ESS A +

2Pn

达到最小。其中ESS(A)表示在子集

A下生成回归方程的残差平方和，P为子集A中变量的个数。

BIC准则：选择子集B，使BIC B =ln ESS B +

Pln(n)n

达到最小。其中ESS(B)表示在子

集B下生成回归方程的残差平方和，P为子集B中变量的个数。

在我的程序里，先生成一个64*6的二进制矩阵，0表示没有选入，1表示选入，该矩阵代表全子集。然后在每一种子集里把AIC和BIC求出来，找出最小的AIC和BIC对应的子集。

3、Stepwise方法

基本思想是逐个引进自变量，每次引入对Y影响最显著的自变量，同时对已选入的变量进行检验，把其中不显著的变量剔除，一直重复此步骤直到方程中全是显著变量且不遗漏该选入的变量。

在我的程序中主要是调用r软件已有的逐步回归函数来实现的。

4、模拟结果：

Lasso算法与AIC、BIC、stepwise方法的比较 估计方法 Lasso AIC BIC stepwise

三、结论：

从错误率的角度看，incorrect.no应该是越接近于0越好，在模拟结果中Lasso算法的错误率是最低的，说明Lasso算法在错误率方面具有较为明显的优势；

从正确率的角度看，correct.no应该是越接近于3越好，在模拟结果中BIC的最高，但同时BIC的incorrect.no也很大，说明用BIC准则估计惩罚力度过大，而Lasso算法的correct比AIC准则、stepwise都高，并且incorrect也低，综上说明Lasso算法具有较好的预测精度，一定程度上优于其它三种方法，可推广使用。

Corret.no 2.585 2.523 2.926 2.51 Incorrect.no 0.201 0.355 1.008 0.348 四、R程序代码：

#为计算四种算法的correct和incorrect分别赋予初值0 corr<-0;incorr<-0 corra<-0;incorra<-0 corrb<-0;incorrb<-0 corrs<-0;incorrs<-0

for(i in 1:1000) {

library(MASS);

A<-rep(1,6); I<-diag(A); B<-rep(0,6);

X<-mvrnorm(200,B,I) #生成随机数自变量X beta<-matrix(c(1,0.8,0.6,0,0,0),6,1) yifu<-rnorm(200,0,4)

Y<-X%*?ta+yifu #生成因变量Y lamida<-seq(0.2,0.5,0.005) #对lam进行取值循环

Bet<-matrix(rep(0,366),6,61) #把大矩阵Bet用来装每个lamida生成的beta

G<-matrix(rep(0,61),1,61) #用G装每个lamida生成的比较式，用于判断大小 for (k in 1:61 ) #对不同的lamida分别进行求相应beta的值 {

lam<-lamida[k]

beta0<-(solve(t(X)%*%X)%*%t(X))%*%Y bi<-matrix(rep(0,6),6,1)

while(t(beta0-bi)%*%(beta0-bi)>1e-6) #当迭代的两次差的模足够小时，停止 {

sgn<-matrix(rep(0,6),6,1) for (k1 in 1:6) {

if (abs(beta0[k1,1]<1e-4)) beta0[k1,1]<-1e-6 #如果beta的绝对值太小，赋予定值 sgn[k1,1]<-1/abs(beta0[k1,1]) }

U2<-(-1)*t(X)%*%(Y-X%*?ta0)+200*lam*diag(c(sgn))%*?ta0; #Q的一次导数U2 U3<-t(X)%*%X+200*lam*diag(c(sgn)); #求Q的二次导 beta1<-beta0-solve(U3)%*%U2 #求出迭代的beta

for(k3 in 1:6) #第三步选择，把小于一个定值的元素设为0 {

if(abs(beta1[k3,1])<1e-1) beta1[k3,1]<-0 }

bi<-beta0; #换元用于下一次同样的迭代 beta0<-beta1;

} #迭代出来是一个lamida生成的beta0

Bet[,k]<-beta0 #将每个lamida生成beta0放入大矩阵Bet