初学线性回归

线性回归是机器学习中的一个非常重要的算法，一般用来进行数据的拟合或者预测，同时还可以用来特征重要性评估。

线性回归一般分为一元线性回归(一个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)