我的R语言小白之梯度上升和逐步回归的结合使用

我们今天的主题通常在用sas拟合逻辑回归模型的时候，我们会使用逐步回归，最优得分统计模型的等方法去拟合模型。而在接触机器学习算法用R和python实践之后，我们会了解到梯度上升算法，和梯度下降算法。其实本质上模型在拟合的时候用的就是最大似然估计来确定逐步回归选出来的一个参数估计，但是这个过程你说看不到，那么现在假设你过程你可以选择，就是你来算这个最大似然估计的过程。甚至，你可以定义这个过程损失函数，那么就要使用最大似然估计。

逐步回归法结合了向前选择法和向后选择法的优点。一开始模型只有截距项，先使用前向选择法选入卡方统计量最大，符合选入模型P值的变量，然后使用后向选择法移除P值最大的变量，即最不显著的变量，不断重复以上过程。所以也可以说逐步回归的每一步都结合了向前选择法和向后选择法。

要学习梯度上升算法和梯度下降算法，就要先了解梯度的概念，要了解梯度就离不开方向导数。学过大学微积分或数学分析的同学都知道，导数代表了一个函数的变化率。但当一个函数包含多个自变量的时候，函数值的变化不仅取决于自变量的变化，还取决于使用哪个自变量。换句话说，函数值同时决定于移动的距离和移动的方向。

然后，梯度其实就是一定最大的方向导数。在自变量只有一个的时候，一点的导数其实是确定的。而到了多个自变量的时候，以一个三维空间为例（如下图的高山），概括为Y为X1，X2的函数，那么在高山上的点上升的方向就不唯一，即方向导数不唯一，那么在某点上山最快的方向就可以描述为该点的梯度。在每爬到一个地方，就不断调整上升最快的方向，最终就可以爬到山顶，成为人生赢家。在算法上就描述为每达到一个移动的步长，就计算该点的梯度，不断使Y值增加，达到最大的Y，最后可以求得最优的X1和X2。

换到梯度下降法，就可以把三维图形想象成一个碗，要想到碗底的话，就应该沿下降最快的方向。数学上就是没一步都求梯度的反方向，最后目标就是求Y的最小值。

说了这么多，那么梯度上升法和下降法对逻辑回归到底有什么用呢？逻辑回归建模有一个目标就是求解最优的系数使似然函数最大化。而下降法可以用来是损失函数最小化。先说似然函数最大化，我们可以令模型的系数为刚才举得例子的x1,x2即自变量，那么我们就可以不断迭代，找到最后的最大的似然函数和最佳的一组系数。系数的梯度上升迭代式可以写为，下面的α就是移动的步长，所乘的就是梯度。

所以，我们可以发现，逐步回归等算法其实优化模型的入模变量，梯度上升法是在选定入模变量之后，求最佳的系数去优化模型。那么，在实践上我们就可以在sas拟合完模型，选定变量后，在用R或者python用梯度上升法去求解最优的系数,但是需要明确一点吗，说是最优那是基于损失函数是一个凸函数，当损失函数不是凸函数的时候，只是找到的是局部最优。L()这个函数是自己定义的一个损失函数组成的一个类似最大似然估计的一个函数。

具体了解下，还是看不懂，可以复习一下导数，偏导数以及方向导数。因为梯度的内容实在有点多，所以还是希望大家对于梯度不了解的，可以上网查询了解。我最初想用这个的时候，我是在想一个问题，就是假设我不用最大似然估计定义的损失函数，假设我想用其他损失函数拟合参数，那我该怎么办，所以才有了今天的分享，可能我说的优点乱，我给出梯度拟合参数的过程，你可能会清晰些：

那么作为R语言小白的我，要出动献出我的梯度上升的代码了，还是参考别人的更改，这里的数据集使用的是你逐步回归选下来的变量。这里这是重新拟合参数，不适用你原来拟合的参数，是不是很作，我也觉得我很作。链接在这：http://blog.csdn.net/yuanhangzhegogo/article/details/40613951。

D<-F[-which(names(F) %in% c('APPL_ID','APPL_STATUS_1'))]

# 为等下产生的样本的矩阵做准备，所以把主键还有因变量删掉

Y=F$APPL_STATUS_1

# 将因变量单独拿出来，等下要进行运算

m<-length(Y)

# 取出y的长度，为的是等下构造截距变量，设为1

#自变量增加一列构造矩阵

x1<-rep(1,m)

# 生成截距变量，设为1

Y<-as.matrix(F$APPL_STATUS_1)

# 生成因变量的矩阵等下可以计算

X<-as.matrix(data.frame(x1,D))

#生成自变量矩阵，等下计算

maxiteration=2000

#设定迭代次数

theta<-matrix(rep(0,14),ncol=1)

# 设置初始的系数

#设定学习速度

alpha=0.0001

pred<-data.frame()

# 生成一个空表

for ( n in c(1:maxiteration)){

#计算梯度

p<-1/(1+exp(-X%*%theta))

#计算通过填入参数之后的预测概率

grad=t(X)%*%(Y-p)

#放入公式计算

a<-theta

# 把前一个的参数矩阵赋给a

theta=theta+alpha*grad

# 计算梯度上升的一个参数

interval<-theta-a

# 计算之间的差值

dd<-data.frame(interval,sum=sum(interval),theta)

# 合并差值，差值的累计，以及对应的参数

pred<-rbind(pred,dd)

# 纵向合并每一次迭代的数据

print(n)

# 打印迭代到哪里，好检查错误以及进度

}

出来的结果看数据集看pred：

最后一列是参数估计，中间是两次梯度相减的累加，可以看到迭代了2000次之后，他的差距已经很小很小的，基本可以断定快到山顶了，你要是觉得这样子差距还是让你不满意，你可以设置迭代次数到3000次。第一列是两个梯度的各个值的相减，这是为了让你看到迭代的过程该变量的权重的变小了还是变大了。当然你也可以更改我的代码，把他改成迭代到两次相减的数小于你设置的数就停止。