所有由郝智恒发布的文章

关于郝智恒

南开数理统计小硕士。主攻计算机试验的设计。

新浪微博文本分析初探v0.1

v0.1版本说明:本文发在主站上之后,站友们经常评论代码跑着有问题。经过和lijian大哥等人进行咨询,自己也摸索了一些之后,发现了之前代码非常多的漏洞。因此,给广大站友带来了困扰。在这里我表示万分的抱歉。最近邮箱中收到让我整理代码的需求越来越多。我也非常想整理下,但是由于工作也非常繁忙,所以很难抽出时间。前两天说5.1期间会整理一下代码发出来。但是事实上因为5.1小长假期间我可能无法上网,导致无力更新代码。所以今晚抽时间对代码进行了简单的修改。对本文也进行了一些调整。目前的状况是,基本上到生成待分析的语料库已经没有问题了。聚类分析的模型可以跑出来,但是最终的图画不出来。我暂时也没能找到原因。所以进一步的调整可能要等到5.1过完以后再抽时间来做了。这篇文章我会负责到底的哈。(20130429)

完整代码如下:weibo_analysis

—————————–分割线————————————————-

自从lijian大哥的Rweibo包问世以来,便成了R爱好者们获取新浪微博数据的最为重要的工具。在该包的中文主页上,作者对如何连接新浪微博的API,获取授权,并以此为基础开发应用的原理讲解的非常清楚。对于我这种连基本的网页开发神马原理都一点也不清楚的菜鸟来说,Rweibo是一种非常趁手的获取微博数据的工具。

有了获取数据的工具,对于中文文本分析来说,最重要的是分词。这里使用的分词算法来自中科院 ictclas算法。依然是沾了lijian大哥Rwordseg的光,直接拿来用了。

有了这两样利器,我们便可以来分析一下新浪微博的数据了。我选取的话题是最近热映的国产喜剧电影《泰囧》,在微博上拿到了998条和“泰囧”有关的微博文本。代码如下(以下代码不能直接执行,请首先阅读链接中Rweibo的关于授权帮助文档):

#关键词搜索并不需要注册API
require(Rweibo)
#registerApp(app_name = "SNA3", "********", "****************")
#roauth <- createOAuth(app_name = "SNA3", access_name = "rweibo")
res <- web.search.content("泰囧", page = 10, sleepmean = 10,
           sleepsd = 1)$Weibo

获取了数据之后,首先迫不及待对微博文本进行分词。代码如下(Rwordseg包可以在语料库中自助加入新词,比如下面的insertWords语句):

require(Rwordseg)
insertWords("泰囧")
n = length(res[, 1])
res = res[res!=" "]
words = unlist(lapply(X = res, FUN = segmentCN))
word = lapply(X = words, FUN = strsplit, " ")
v = table(unlist(word))
v = sort(v, decreasing = T)
v[1:100]
head(v)
d = data.frame(word = names(v), freq = v)

继续阅读新浪微博文本分析初探v0.1

计算机试验简介

很早就想为COS写一篇关于计算机试验的东西。可是始终也未敢动笔,觉得自己才疏学浅,生怕写得偏颇。但是另外一方面,又觉得这是COS上一块空白的话题,没有人提及过。今天写这篇文章,主要是为抛砖引玉,另外丰富一下主站文章的题材。

在这篇文章中,我想大致介绍一下计算机试验的设计以及建模,另外会有一些R中专门做计算机试验的包的相关介绍。

目前,比较流行的计算机试验设计与建模的教材有两本(似乎也只有这两本)。一本是Thomas J. Santner, Brian J.Williams以及William I.Notz合著的,Springer2003年出版的《The Design and Analysis of Computer Experiment》。另外一本是Kai-Tai Fang ,Runze Li以及Agus Sudjianto合著的,2005年由Chapman&Hall出版的《Design and Modeling for Computer Experiment》。所谓的计算机试验(computer experiment),是相对于传统的实体试验(physical experiment)而言的。我们都知道,实体试验的数据(比如农业试验,生物试验等得到的数据)总是会受到随机误差的影响,因此在实体试验设计中,引入了“重复”的原则,目的就是要通过重复的试验,减小随机误差对于分析的影响。(关于实体试验的设计以及建模,国内外已经有很多教材讨论这些,我相信大多数统计专业的同学所学习的试验设计入门课程都是在学习实体试验的设计以及数据分析。)与此相对的,计算机试验得到的数据并不受到随机误差的干扰。因为一段固定的代码,在计算机上无论运行多少次,得到的结果都是一样的。

计算机试验,通常比实体试验包含了更多的变量,用于研究特别复杂的系统,比如航天探测器。对于这些特别复杂的系统,往往需要用一个更加简单的拟模型(meta model)来逼近。计算机试验的目的之一,就是要通过在计算机上,利用代码模拟某个复杂的设备或者某个过程,得到了试验的数据之后,通过分析建模,从而得到一个相对更简单的模型,被称为拟模型,拟模型是非常有用的。

由于计算机试验具有复杂性,以及输出结果的确定性,做计算机试验收集数据的时候,我们需要一些特殊的设计。

试验设计,总是和模型紧密相联系的。事实上,试验设计是一种收集数据的策略。在实体试验中,当我们需要比较某种因子不同水平之间的差别,我们通常会选用因析设计。因析设计事实上就是和因子模型相对应而产生的收集数据的策略。另外,所谓的最优设计,也是相对应于模型的最优化准则而产生的收集数据的策略。而在计算机试验中,常用的是模型是所谓的全局均值模型(overall mean model)。

不失一般性,我们考虑试验区域是一个s维的单位立方$C^s=[0,1]^s$。当试验次数给定了,比如说$n$次,我们就需要考虑如何找到一个好的试验$D_n=\{x_1,\cdots,x_n\}$,使得$f(x)-g(x)$对于试验区域上的所有点,都尽可能的小。其中$f(x)$是真实的模型,而$g(x)$是拟模型。我们考虑到,对于所有的点,都要求上述差分达到最小是很困难的。因此,很多学者都将该问题进行了简化,转而去寻找全局均值$E(y)=\int_{C^s}f(x)dx$的最佳的估计值。而通常,都用样本均值,也就是平均数$\bar{y}(D_n)=\frac{1}{n}\sum_{i=1}^nf(x_i)$来估计总体均值。因此,我们的问题也就转化为了如何找到一个试验,使得这个试验得到的数据做样本均值来估计总体均值是最好的。这也就是计算机试验设计的动因。

从统计学的角度来看,如果样本点$x_1,\cdots,x_n$是在$C^s$上均匀分布上独立抽取的,那么样本均值是无偏的,并且方差为$\frac{var(f(x))}{n}$。但是,通常来说,这个方差太大了。因此很多学者提出了不同的抽样方法来降低这个样本均值的方差。1979年,McKay,Beckman和Conover基于分层抽样提出了一种新的抽样方法,这种方法使得随机选择的样本点并不独立,并且边际分布相同,因此

$Var(\bar{y}(D_n))=\frac{1}{n}Var(f(x))+\frac{n-1}{n}Cov(f(x_1),f(x_2))$

当$f(x_1)$与$f(x_2)$负相关时,方差就减小了。这种方法被人不断的推广,创新。成为了现在计算机试验设计中一块非常重要的内容——拉丁超立方体设计(Latin-Hypercube design)。在后来的不断发展中,拉丁超立方和正交表相结合,产生了很多正交拉丁超立方设计。此外,当引入一些准则时,比如熵准则,最大最小距离准则等,也产生了对应该准则的最优的设计,包括最大最小拉丁超立方等。这些设计在计算机试验中应用非常广泛,最近也有人将这种设计用在变量选择的过程中,通过模拟,大大提高了变量选择的精确度。

与拉丁超立方设计相对应的另外一种设计,是由我国科学家方开泰,王元在上世纪八十年代提出的均匀设计。均匀设计这一想法的提出,起源于伪蒙特卡洛方法中的Koksma-Hlawaka不等式。$|E(y)-\bar{y}(D_n)|\leq V(f)D(D_n)$。其中$D(D_n)$是设计$D_n$的星偏差(star-discrepancy),而$V(f)$是方程$f$的总体变动情况。其中星偏差是一种均匀性的度量,当试验点在试验区域分布越均匀时,星偏差的值就越小,因此,样本均值与总体均值之差的绝对值的上界也就越小,相当于提高了估计的精度。一个试验区域上具有最小星偏差的设计,就是所谓的均匀设计。均匀设计由我国科学家自主提出,并且在国际上得到了非常好的反响以及应用。之后,有很多学者,对于均匀设计做出了改进,使得均匀设计也不断地发展壮大。

下面简单介绍一下对于计算机试验得到的数据建模的常用方法。事实上,统计学习中的常用的建模方法都可以用来对计算机试验数据进行建模。包括样条回归,局部回归,神经网络等,都可以对数据进行建模。这里要要重点介绍的是一种非常常用的模型,高斯-克里金模型(Gaussian-Kriging models)。克里金是南非地质学家,他在其硕士论文中分析矿产数据的时候,提出了克里金模型,后来由几位统计学家发展了他的模型,得到了目前在计算机试验中非常常用的高斯-克里金模型。

$y(x)=\sum_{j=0}^{L}\beta_jB_j(x)+z(x)$

其中$B_j(x)$是选定的基函数,$z(x)$是随机误差,来自于高斯过程。该高斯过程零均值,方差为$\sigma^2$,协方差函数为$r(\theta,s,t)=corr(z(s),z(t))$。通常情况下$r$是预先选定的。可以证明,高斯克里金模型的预测值是BLUP。

关于计算机试验数据的建模,方法很多,但是建模方法仅仅是建模方法,无论是计算机试验还是实体试验,建模方法往往通用,因此,这里不做展开论述了。

下面简单介绍一下R语言中,做计算机试验的一些包。

lhs包,主要用来生成拉丁超立方体设计。DiceDesign包可以生成更多的一些空间填充设计。DiceKriging提供了对计算机试验数据构建高斯-克里金模型的功能,DiceView则提供了多维模型可视化的手段。以上的包在task view中都有所介绍。此外,值得一提的是VizCompX包,其对计算机试验的数据拟合高斯模型,并且提供可视化的手段。感兴趣的童鞋可以参看这几个包的帮助文档,帮助文档都不长,容易上手。

最后,提一点我个人的小感想。我认为,计算机试验中,更应该引人注意的话题应该是试验的设计。一个好的设计,不光能够提高试验的精度,更能体现一种美感。生成一个好的设计,大致上有两种路可以走,一种是预先设定一个准则,然后利用最优化的算法,进行搜索。常用的算法包括一些全局的最优算法,比如门限接受算法,模拟退火算法等。另外一种方法,则是利用数学中其他学科,比如伽罗瓦理论,编码理论等,通过这些学科的理论特点,精巧地构造一个设计。这两种构造设计的方法相互补充,相辅相成,提供了非常广大的研究空间。另外,这些为计算机试验而提出的设计,也可以用在很多其他的统计学的课题中,比如之前已经提到过的变量选择方法精确度的提升方面。总之,试验设计是研究如何有效收集数据的,无论数据是来自于实体试验,还是计算机试验,一个设计了的试验,往往都可以提高效率,提高精度。这也是试验设计这门课程的意义所在。越南的Nyuan开发了一款试验设计的软件Gendex,通过输入试验的行数,列数,选择你所需要的设计种类,便可以生成多种不同的设计阵。这款软件,是目前可以生成设计阵种类最多的,感兴趣的童鞋,不妨搜来玩一玩。

分组最小角回归算法(group LARS)

继续前两篇博文中对于最小角回归(LARS)和lasso的介绍。在这篇文章中,我打算介绍一下分组最小角回归算法(Group LARS)。本文的主要观点均来自Ming Yuan和Yi Lin二人2006合作发表在JRSSB上的论文Model selection and estimation in regression with grouped variables.

首先,我想说明一下,为何要引入分组变量(grouped variable)的概念。举一个简单的例子,在可加的多项式模型中,每一项都是多项式。这个多项式有可能可以通过最初的变量的线性组合来表达。在进行这种类型的回归中,挑选重要的变量其实质是挑选重要的因子(factor),而因子则是最初的那些变量的线性组合。分组变量的回归问题,实际上就是我们一般所说的回归问题的推广。如果我们把每一个单独的变量都看成一个因子,那么这种情况下的回归就是一般意义下的回归。下面用公式更加直白的说明这个问题:

$Y=\sum_{j=1}^JX_j\beta_j+e$

其中$Y$是个$n$维向量,$e\~N_n(0,\sigma^2I)$.$X_j$是$n\times p_j$矩阵,代表的是第j个因子(factor,是变量variables的线性组合)。$\beta_j$是$p_j$维的系数向量。依然假定$Y$是中心化的,$X_j$是中心化并且正交化的($X_j’X_j=I$)。这个就是分组变量的回归模型。

在小弟之前的文章中介绍过最小角回归(LARS)算法,对于变量选择来说,这种算法是比较易于计算机实现的,而且比传统的向前逐步回归(forward selection)更加谨慎,但是没有逐段回归(stagewise)那么谨小慎微。对于分组变量的回归模型选择问题来说,直接套用LARS算法可能会挑选出一些不必要的因子进入回归模型中。因此在文章中,作者提出了分组最小角回归(group LARS).这种算法对最小角回归进行了推广,下面简单介绍一下这种算法。

推广是循序渐进的,因此,我们先来看看$p_j$均相等的情况(注意这里$p_j$未必等于1)。在传统的LARS算法中,我们是要保证入选变量和当前残差的角度都相同。如今,我们就需要保证入选的因子们(再次强调,这里是因子,是原始变量的线性组合)与当前残差的角度相同。直观地来定义一个角度$\theta(r,X_j)$,这个角度就是向量$r$与$X_j$中的变量所构成的空间的夹角,也就是$r$与它在$X_j$中的变量所构成空间的投影的夹角。(再废话两句,这个夹角越小,就意味着$r$与$X_j$中的变量构成的空间中的相关系数越大。)当$r$是当前残差的时候,我们就可以确定求解路径(solution path)为当前残差在$X_j$中的变量所构成的空间上的投影的方向。

上面说的可能还不够清楚,下面我分步骤具体介绍每一步的做法。

  1. 最开始时,回归变量集依然是空集,此时的残差就是$Y$,我们寻找和$Y$夹角最小的$X_j$.记它为$X_{j1}$.(这里要插一句,就是如何度量这个角度。用这个角的余弦值的平方即可:$cos^2(\theta(r,X_j))=\frac{\|X_j’r\|^2}{\|r\|^2}$,角度越小,这个值越大。)
  2. 确定求解路径为$Y$在$X_{j1}$中的变量所构成的空间上的投影的方向。然后在这个路径上前进,直到出现第二个因子$X_{j2}$,使得当前残差$r$与这两个因子的夹角同样小。也就是$\|X_{j1}’r\|^2=\|X_{j2}’r\|^2$。
  3. 计算新的求解路径,这个路径其实就是当前残差在$X_{j1}$和$X_{j2}$的变量所构成的空间中投影的方向。然后再沿着这条路径前进,直到第三个因子出现,使得当前残差和这三个因子的夹角都是同样小。

然后继续按照上面的步骤,计算出求解路径,然后继续前进。直到找到第四个满足条件的因子。

这些步骤和原始的LARS的进程是同样的。只不过在确定求解方向时有了变化。

上面介绍的都是$p_j$相同的情况,当$p_j$不同的时候,只要对上面的过程做一个小的改动即可。用$\frac{\|X_j’r\|^2}{p_j}$代替$\|X_j’r\|^2$,其他过程都不变。这种对于LARS算法的推广是比较直观的。而且M.Yuan和Y.Lin在文章中还指出,这种算法相对于正交化$X_j$的过程是稳定的,通过Gram-Schmidt方法进行正交化本来是依赖于选择的参数的,但是分组LARS确定的求解方向并不依赖该参数。

修正的LARS算法和lasso

在小弟的上一篇文章中,简单的介绍了LARS算法是怎么回事。主要参考的是Efron等人的经典文章least angle regression。在这篇文章中,还提到了一些有趣的看法,比如如何用LARS算法来求解lasso estimate和forward stagewise estimate。这种看法将我对于模型选择的认识提升了一个层次。在这个更高的层次下看回归的变量选择过程,似乎能有一些更加创新的想法。

lasso estimate的提出是Tibshirani在1996年JRSSB上的一篇文章Regression shrinkage and selection via lasso。所谓lasso,其全称是least absolute shrinkage and selection operator。其想法可以用如下的最优化问题来表述:

在限制了$\sum_{j=1}^J|\hat{\beta_j}|\leq t$的情况下,求使得残差平和$\|y-X\hat{\beta}\|^2$达到最小的回归系数的估值。

我们熟悉如何求解限制条件为等号时,回归方程的求解。也就是用lagrange乘子法求解。但是对于这种,限制条件是不等号的情况,该如何求解,则有两种想法。第一种,也是我比较倾向于的方法,是利用计算机程序,对$t$从$0$开始,不断慢慢增加它的值,然后对每个$t$,求限制条件为等号时候的回归系数的估计,从而可以以$t$的值为横轴,作出一系列的回归系数向量的估计值,这一系列的回归系数的估计值就是lasso estimation。另外一种想法,是借助与最优化问题中的KKT条件,用某个黑箱式的算法,求解。(本人对于最优化方面的东西实在是不很熟悉,故不在此弄斧,只求抛砖引玉,能有高手给出这种想法的具体介绍。)

lasso estimate具有shrinkage和selection两种功能,shrinkage这个不用多讲,本科期间学过回归分析的同学应该都知道岭估计会有shrinkage的功效,lasso也同样。关于selection功能,Tibshirani提出,当$t$值小到一定程度的时候,lasso estimate会使得某些回归系数的估值是$0$,这确实是起到了变量选择的作用。当$t$不断增大时,选入回归模型的变量会逐渐增多,当$t$增大到某个值时,所有变量都入选了回归模型,这个时候得到的回归模型的系数是通常意义下的最小二乘估计。从这个角度上来看,lasso也可以看做是一种逐步回归的过程。

在我的上一篇文章中,提到了Efron对于逐步回归的一种看法,就是在某个标准之下(比如LARS的标准就是要保证当前残差和已入选变量之间的相关系数相等,也就是当前残差在已入选变量的构成空间中的投影,是那些变量的角平分线)选择一条solution path,在这个solution path上proceed,不断吸收新的变量进入,然后调整solution path 继续proceed。那么对于求解lasso的算法,也有一个相应的对应。Efron提出了一种修正的LARS算法,可以用修正的LARS算法来求解所有的lasso estimates。下面我介绍一下这种修正的LARS算法。

首先假设我们已经完成了几步LARS steps。这时候,我们已经有了一个回归变量集,我们记这个回归变量集为$X_{\mathscr{A}}$。这个集合就对应着一个对于$Y$的估计,我们记为$\hat{\mu}_{\mathscr{A}}$。这个估值对应着一个lasso方法对于响应的估值(这里我认为LARS估值和lasso估值应该是一样的),lasso的估值,对应着回归系数的lasso估值,回归系数向量的lasso估值我们记为$\hat{\beta}$。

为了继续进行下一步,我们先给出一个向量的表达式,然后再解释一下它

$w_{A}=(1_{A}'(X_{A}’X_{A})^{-1}1_{A})^{-\frac{1}{2}}(X_{A}’X_{A})^{-1}1_{A}$.

$X_{A}w_{A}$就是LARS算法的在当前回归变量集下的solution path。那么我们可以把$w_{A}$作为$\beta$的proceed的path。Efron定义了一个向量$\hat{d}$,这个向量的元素是$s_jw_j$,其中$s_j$是入选变量$x_j$与当前残差的相关系数的符号,也是$\hat{\beta_j}$的符号。对于没有入选的变量,他们对应在$\hat{d}$中的元素为0。也就是对应着$\mu(r)=X\beta(r)$,我们有

$\beta_j(r)=\hat{\beta_j}+r\hat{d_j}$

将LARS的solution path对应到lasso estimate的path上,这种对应的想法非常值得借鉴。

很显然,$\beta_j(r)$会在$r_j=-\hat{\beta_j}/\hat{d_j}$处变号。那么对于我们已经有的lasso estimate$\beta(r)$,它中的元素会在最小的的那个大于$0$的$r_j$处变号。我们记之为$\bar{r}$。如果没有$r_j$大于$0$,那么$\bar{r}$就记为无穷大。

对于LARS本身而言,在已经有了如今的回归变量集和当前残差的基础上,我们就会有条solution path,在这个solution path上proceed的最大步记为$\hat{r}$.通过比较$\hat{r}$和$\bar{r}$就会有进一步的想法。Efron的文章证明了如果$\bar{r}$小于$\hat{r}$,则对应于LARS估计的那个$\beta_j(r)$不会成为一个lasso estimation。(这个是因为当前残差和对应变量的相关系数的符号一定是和该变量的系数符号一致才行)。在这种情况下,我们就不能继续在LARS的solution path上继续前进了,为了利用LARS算法求得lasso estimate,Efron提出把$\bar{r}$所对应的那个$r_j$所对应的$x_j$从回归变量中去掉。去掉之后再计算当前残差和当前这些变量集之间的相关系数,从而确定一条新的solution path,继续进行LARS step。这样进行下去,可以通过LARS算法得到所有的lasso estimate。

这个对于LARS的lasso修正算法,被Efron称作“one at a time”条件,也就是每一步都要增加或删掉一个变量。下图显示了用修正了的LARS算法求lasso estimate的过程。

这个图是Efron等人的文章中,对于一个实际数据进行回归得到的。该数据一共有10个变量。图的横轴,是所有回归系数估值的绝对值之和,这个值从$0$增加。左侧的纵轴,是回归系数的估值,右侧纵轴是这些回归系数对应的变量的下标。这个图中,我们可以看到每一个回归系数的path。可以看到第七个变量对应的回归系数在横轴快到3000的时候变为了0,说明到这一步时,该变量被删除掉,之后又被重新添加到了回归变量集中。

下面通过一个简单的模拟,对lars和lasso以及forward stagewise做一个简单的实现。其实在R中已经有了一个名为lars的包,可以实现上述三种回归。

首先,我要模拟的方程为

$y={x^3}_1+{x^2}_1+x_1+\frac{1}{3}{x^3}_2-{x^2}_2+\frac{2}{3}x_2+e$

其中$x_1$和$x_2$是服从二维联合正态分布,均值为零向量,$cov(x_1,x_2)=0.5$,$var(x_1)=var(x_2)=1$,$e$服从$N(0,9)$。我取了50次观测,然后分别通过lasso,lars,以及forward stagewise三种算法进行了回归,其变量的回归路径如下图。

简单的代码我直接贴在本文的最后。从这三个算法的图中,我们并看不出有特别的区别,只能看出一些细小的差别。至于要判断哪种算法更好,则应该因问题而异。也不是本文能够论述的问题了。

对于LARS算法的修正,还可以应用到计算forward stagewise的estimate中,在Efron的文章中也有介绍。他的这种看法,好似凌驾在整个回归变量选择过程之上,从一个更高的角度观察之,给出一种更为一般性的视角。这也就是大牛和一般人之间的差别。读Efron的文章,总有一种让人想要膜拜的冲动。对于模型选择方面的东西,值得挖掘的还很多。Tibshirani在最新的一篇综述性的文章中,给出了lasso的诞生到现今发展的一系列流程。感兴趣的读者,可以去看看这篇文章,在cos论坛上有。链接如下:

http://cos.name/cn/topic/104104

用lars算法做模拟的代码:

利用lars模拟

 

LARS算法简介

最近临时抱佛脚,为了讨论班报告Group Regression方面的文章,研究了Efron等人于2004年发表在Annals of Statistics里一篇被讨论的文章LEAST ANGLE REGRESSION。这篇文章很长,有45页。加上后面一些模型方面大牛的讨论的文章,一共有93页。对于这种超长论文,我向来敬畏。后来因为要报告的文章里很多东西都看不懂,才回过头来研读这篇基石性的文章。

所谓大牛,就是他能提出一种别人从来没有提出过的想法。大牛们看待问题的角度和常人不同。比如在回归中常用的逐步回归法。我们小辈们只知道向前回归,向后回归还有二者结合的一些最基本的想法。比如向前回归,就是先选择和响应最相关的变量,进行最小二乘回归。然后在这个模型的基础上,再选择和此时残差相关度最高的(也就是相关度次高)的变量,加入模型重新最小二乘回归。之后再如法继续,直到在某些度量模型的最优性准则之下达到最优,从而选取一个最优的变量子集进行回归分析,得到的模型是相比原模型更加简便,更易于解释的。这种方法,牺牲了模型准确性(预测有偏),但是提高了模型的精确度(方差变小)。大多数本科生对逐步回归的理解也就如此了。Efron看待这个问题时,比起常人更高了一个层次。他首先指出,逐步向前回归,有可能在第二步挑选变量的时候去掉和X1相关的,但是也很重要的解释变量。这是因为它每次找到变量,前进的步伐都太大了,侵略性太强。

因此在这个基础上,Efron提出了Forward stagewise。也就是先找出和响应最相关的一个变量,找到第一个变量后不急于做最小二乘回归,而是在变量的solution path上一点一点的前进(所谓solution path是指一个方向,逐步回归是在这个方向上进行),每前进一点,都要计算一下当前的残差和原有的所有变量的相关系数,找出绝对值最大的相关系数对应的变量。我们可以想像,刚开始,前进的步伐很小,相关系数绝对值最大的对应的变量一定还是第一步选入的变量。但是随着前进的进程不断向前,这个相关系数的绝对值是在慢慢减小的,直到找到另外一个变量X2,它和当前前残差的相关系数和第一个入选变量X1的相关系数绝对值相同,并列第一。此时把X2也加入回归模型中,此时回归模型在X1上的系数已经确定了,如果在X1的solution path上继续前进,则得到的与当前残差相关系数最大的变量一定是X2,所以不再前进,而是改为在X2的solution path上前进,直到找到第三个变量X3,使得X3的与当前残差的相关系数绝对值最大。这样一步一步进行下去。每一步都是很多小步组成。直到某个模型判定准则生效,停止这个步骤。在每一个solution path上的计算都是线性的。总体的solution path是分段线性的。这种算法是一种自动进行模型构建的方法。它和传统的Forward selection在本质上是一样的,都是选择一个变量,然后选择一个继续进行的solution path,在该方向上前进。这两种方法的solution path的选择方法是一样的,唯一的区别就是前进的步伐不一样,Forward selection的前进步伐很大,一次到头,而stagewise则是一小步一小步前进。这样比Forward selection要谨慎一些,会免于漏掉一些重要的变量。

从这个视角来看,我们可以选择另外一种solution path。Efron等人在这篇文章中,就提出了一种新的solution path。在已经入选的变量中,寻找一个新的路径,使得在这个路径上前进时,当前残差与已入选变量的相关系数都是相同的。直到找出新的与当前残差相关系数最大的变量。从几何上来看,当前残差在那些已选入回归集的变量们所构成的空间中的投影,是这些变量的角平分线。下面我简单的描述一下这个算法:

  • 第一步,我们初始的估计模型为0,那么当前的残差就是Y,我们找出X’Y中绝对值最大的那个对应的变量,记为X1,把它加入回归模型。这一步中X’Y是当前残差和所有变量的相关系数向量。(注意这里Y都已经中心化,X中心标准化过了)。
  • 第二步,在已选的变量的solution path上前进,solution path就是s1*X1,s1是X1与当前残差的相关系数的符号。在这个path上前进,直到另外一个变量出现,使得X1与当前残差的相关系数与它和当前残差的相关系数相同。记这个变量为X2,把它加入回归模型中。
  • 第三步,找到新的solution path。Efron在文章中提出了一种找出满足LARS条件的solution path的解法。solution path需要使得已选入模型变量和当前残差的相关系数均相等。因此这样的路径选择它的方向很显然就是$X_k(X_k’X_k)^{-1}1$的指向(因为$X_k'(X_k(X_k’X_k)^{-1})1$的元素都相同,保证了LARS的要求,当然这里或许会有一些其他的解,也能满足LARS的要求,有没有达人能想到或许证明这个解是唯一的)。只要再标准化这个向量,我们便就找到了solution path的方向。在这个方向上前进,直到下一个满足与当前残差相关系数绝对值最大的变量出现。如此继续下去。

LARS算法,保证了所有入选回归模型的变量在solution path上前进的时候,与当前残差的相关系数都是一样的。这一点,比起Forward stagewise要捷径一些,走得更快一些。

LARS算法已经在SAS和R中实现了。作为回归模型选择的一种重要的算法,LARS相比起传统的Forward selection和Forward stagewise,既不那么富于侵略性,又比较走捷径。LARS算法在lasso 估计的求解中也有非常好的应用。在Efron等人的同篇论文中有详细的讨论。关于lasso和它的LARS算法,笔者将在今后的文章中介绍。