线性回归是机器学习中的一个非常重要的算法,一般用来进行数据的拟合或者预测,同时还可以用来特征重要性评估。
线性回归一般分为一元线性回归(一个x,一个y)、多元线性回归(多个x,一个y):
例如 :y=a*x+b,这就是一元函数
y=a*x1+b*x2+c ,这就是多元线性回归
在这里我只讲述一元的线性回归,多元的类型,我会在后面提及
首先,我们这里有一堆数据
所谓一元线性回归要做的事情,实际上是去找一条最优直线,能够以最小的误差去拟合这些点。
如上图所示,横坐标表示x,纵坐标表示y。我们要找的就是图中的这条直线。我们要去找到这条直线,大家可以想象,我们肯定希望找到的那条线,距离每个点都很近,最好所有的点上都在这条线上,但是一条直线去拟合所有的点都在这条直线上肯定不现实,所以我们希望这些点尽量离这条直线近一点。即去找每个点和直线的距离 最小的那条线,一般我们有几个确定参数的方法:
- 使得误差的最大绝对值达到最小
- 使误差绝对值的平均值达到最小
- 使误差的平方和达到最小
第一种方法中,所确定的参数使得最大误差达到最小,因而称为最大最小问题,也称为最优一致逼近问题。这类问题在求解上存在一定的困难
第二种方法中,由于函数关于参数不可微,因此不能用一般的极值方法求最优解
第三中方法中,所确定的参数使得均方差最小,这种逼近称为最小二乘逼近,此时成y=mx+b为最小二乘拟合直线,并且由于函数关于参数连续可微,所以计算上易于处理。在这里,我们采用第三种方法,那么误差(代价函数)可以表示为,这里i表示第i个数据,N表示总的样本个数。一般我们还会把Loss求和平均,来当作最终的损失,即:
怎么求极值呢?
令每个变量的偏导数为零,求方程组的解呗,这个是很基础的高数问题了。
我们可以得到下面的方程组
然后就是巴拉巴拉巴拉把m和b求出来,这样就得到我们要的线性方程了。
梯度下降法
没有梯度下降就没有现在的深度学习,这是一个神奇的算法。
最小二乘法可以一步到位,直接算出m和b,但他是有前提的,具体我有点记不清了,好像是需要满秩什么的。梯度下降法和最小二乘不一样,它通过一步一步的迭代,慢慢的去靠近到那条最优直线。
最小二乘法里面我们提到了两个偏导数,分别为
我们要去找Loss这个方程的最小值,最小值怎么求?按数学的求法就是最小二乘法呗,但是大家可以直观的想一下,很多地方都会用一个碗来形容,那我也找个碗来解释吧。
大家把这个Loss函数想象成这个碗,而我们要求的最小值就是碗底。假设我们现在不能用最小二乘法求极小值,但是我们的计算机的计算能量很强,我们可以用计算量换结果,不管我们位于这个碗的什么位置,只要我们想去碗底,就要往下走。
往下走????????
这个下不就是往梯度方向走吗,那我们沿着梯度一点一点滑下去呗,反正计算机不嫌累。梯度不就是上面那两个公式呗。现在梯度有了,那每次滑多远呢,一滑划过头了不久白算半天了吗,所以还得定义步长,用来表示每次滑多长。这样我们就能每次向下走一点点,再定义一个迭代值用来表示滑多少次,这样我们就能慢慢的一点点的靠近最小值了,不出意外还是能距离最优值很近的。
顺便把上面这个梯度下降法实现下
每次向下滑要慢慢滑,就是要个步长,我们定义为learning_rate,往往很小的一个值。
向下滑动的次数,就是迭代的次数,我定义为num_iter,相对learning_rate往往很大。
定义好这两个,我们就可以一边求梯度,一边向下滑了。就是去更新m和b。
我这里用了加号,很多人会误以为梯度下降就要减,但是其实梯度本身是有方向的,所以这里直接加就可以。
整体代码如下:
import numpy import pylab import math #生成数据集 def createData(): dataSet=numpy.loadtxt('DataSet/Data.txt') #读取Data.txt文件 dataSet=dataSet.astype(dtype=float) #转化成float类型 return dataSet def line_regression(dataSet): #目标函数设为 y=ax+b learn_rate=0.001 #学习速率 init_a=0 #a初始化为0 init_b=0 #b初始化为0 iter_count=1000 #迭代次数 a,b=optimization(dataSet,init_a,init_b,learn_rate,iter_count) #寻找最优的参数a、b,使得拟合效果最佳 print(a,b) plot_data(dataSet, a, b) def plot_data(data,a,b): #plottting x = data[:,0] y = data[:,1] y_predict = a*x+b pylab.plot(x,y,'o') #pylab.plot(x,y_predict,'k-') pylab.show() def optimization(dataSet,init_a,init_b,learn_rate,iter_count): a=init_a b=init_b for i in numpy.arange(0,iter_count): #迭代次数 a,b=compute_gradient(a,b,dataSet,learn_rate) if(i%100==0): print('iter {0}:error={1}'.format(i, compute_error(a, b, dataSet))) return (a,b) def compute_error(a,b,dataSet): sum=0 for i in numpy.arange(len(dataSet)): x=dataSet[i,0] y=dataSet[i,1] y_predict=a*x+b sum+=math.pow(y-y_predict,2) return sum/len(dataSet) def compute_gradient(a,b,dataSet,learn_rata): a_gradient = 0 #a的偏导数 b_gradient = 0 #b的偏导数 N = float(len(dataSet)) # Two ways to implement this # first way for i in range(0,len(dataSet)): #循环每一对数据,求a,b的偏导数 x = dataSet[i,0] y = dataSet[i,1] #computing partial derivations of our error function #b_gradient = -(2/N)*sum((y-(m*x+b))^2) #m_gradient = -(2/N)*sum(x*(y-(m*x+b))^2) a_gradient += -(2 / N) * x * (y - ((a * x) + b)) b_gradient += -(2/N)*(y-((a*x)+b)) # Vectorization implementation # x = dataSet[:, 0] # y = dataSet[:, 1] # b_gradient = -(2 / N) * (y - m_current * x - b_current) # b_gradient = np.sum(b_gradient, axis=0) # m_gradient = -(2 / N) * x * (y - m_current * x - b_current) # m_gradient = np.sum(m_gradient, axis=0) # update our b and m values using out partial derivations #print(a_gradient,b_gradient) new_a = a - (learn_rata * a_gradient) #如果a_gradient>0,说明在a点是递增的,所以按照梯度下降原理,我们要往相反方向走,所以要减 #如果a_gradient<0,说明在a点是递减的,所以按照梯度下降原理,要按照这个方向走下去,所以减去一个 #负值,相当于加上一个正值 new_b = b - (learn_rata * b_gradient) #原理同上 return [new_a, new_b] if __name__=="__main__": dataSet=createData() line_regression(dataSet)