zoukankan      html  css  js  c++  java
  • 大叔学ML第二:线性回归

    基本形式

    线性回归非常直观简洁,是一种常用的回归模型,大叔总结如下:

    设有样本(X)形如:

    [egin{pmatrix} x_1^{(1)} & x_2^{(1)} & cdots &x_n^{(1)}\ x_1^{(2)} & x_2^{(2)} & cdots & x_n^{(2)}\ vdots & vdots & vdots & vdots\ x_1^{(m)} & x_2^{(m)} & cdots & x_n^{(m)}\ end{pmatrix} ]

    对应的标记(vec{y})形如:

    [egin{pmatrix} y^{(1)} \ y^{(2)} \ vdots \ y^{(m)} \ end{pmatrix}]

    其中,矩阵(X)的每一行表示一个样本,一共有m个样本;每列表示样本的一个属性,共有n个属性。设假设函数

    [h(x_1,x_2 dots x_n)= heta_0 + heta_1 x_1 + heta_2 x_2 + dots + heta_n x_n ag{1} ]

    (x_0=1),则(1)式重新写为

    [h(x_1,x_2 dots x_n)= heta_0x_0 + heta_1 x_1 + heta_2 x_2 + dots + heta_n x_n ag{2} ]

    定义代价函数(均方误差)

    [j( heta_0, heta_1dots heta_n)=frac{1}{2m}sum_{k=1}^m (h(x_1^{(k)},x_2^{(k)} dots x_n^{(k)}) - y^{(k)})^2 ]

    即:

    [j( heta_0, heta_1dots heta_n)=frac{1}{2m}sum_{k=1}^m ( heta_0x_0^{(k)} + heta_1 x_1^{(k)} + heta_2 x_2^{(k)} + dots + heta_n x_n^{(k)} - y^{(k)})^2 ag{3} ]

    这里的分母乘以2并没有意义,只是为了求导后正好约掉。另外,其实求绝对值之和更直观,但是计算不方便,求平方后再求和效果是一样的,而且计算非常容易。我们的目标是根据样本数据求出使得代价函数取值最小的参数(vec heta),均方误差越小,说明以(vec heta)为参数的线性函数拟合样本的能力越强

    求解参数(vec heta)

    梯度下降法

    关于梯度下降法可参考 大叔学ML第一:梯度下降

    由于代价函数是一个凸函数,可以用梯度下降法找到最小值。由于用到梯度,首先对( heta_0)( heta_1)( heta_2)直到( heta_n)求偏导:

    • (frac{partial}{partial heta_0}j( heta_0, heta_1dots heta_n) = frac{1}{m}sum_{k=1}^m( heta_0x_0^{(k)} + heta_1x_1^{(k)} + dots+ heta_nx_n^{(k)} - y^{(k)})x_0^{(k)})
    • (frac{partial}{partial heta_1}j( heta_0, heta_1dots heta_n) = frac{1}{m}sum_{k=1}^m( heta_0x_0^{(k)} + heta_1x_1^{(k)} + dots+ heta_nx_n^{(k)}- y^{(k)})x_1^{(k)})
    • (dots)
    • (frac{partial}{partial heta_n}j( heta_0, heta_1dots heta_n) = frac{1}{m}sum_{k=1}^m( heta_0x_0^{(k)} + heta_1x_1^{(k)} + dots+ heta_nx_n^{(k)}- y^{(k)})x_n^{(k)})

    可归纳为:(frac{partial}{partial heta_n}j( heta_0, heta_1dots heta_n) = frac{1}{m}sum_{k=1}^m( heta_0x_0^{(k)} + heta_1x_1^{(k)} + dots+ heta_nx_n^{(k)}- y^{(k)})x_n^{(k)} ag{4})

    万事俱备,现在可以编程了。创建一组测试数据,每组数据包括3个属性,我们来编码拟合出一个线性函数:

    import numpy as np
    
    def gradient(X, Y, m, theta):
        ''' 求theta位置的梯度.
    
        Args:
            X: 样本
            Y: 样本标记
            m: 样本数
            theta: 欲求梯度的位置
        
        Returns:
            gi: theta处函数的梯度值
        '''
        theta_size = np.size(theta)
        g = np.zeros(theta_size)
    
        for i in range(theta_size):   
            gi = 0 #第i个theta分量对应的偏导     
            for j in range(m):
                gi += ((np.dot(X[j], theta) - Y[j]) * X[j, i])
            gi = gi / m 
            g[i] = gi
    
        return g
    
    def gradient_descent(X, Y, step = 0.02, threshold = 0.01):
        ''' 梯度下降法求使代价函数最小的 theta
    
        Args:
            X: 样本
            Y: 样本标记
            step:步长
            threshold:梯度模长阈值,低于此值时停止迭代
        Returns:
            theta: 使代价函数取最小值的theta
        '''
        theta = np.random.rand(4)    
        grad = gradient(X, Y, np.size(X, 0), theta)
        norm = np.linalg.norm(grad)
        
        while(norm > threshold):
            theta -= step * grad
            grad = gradient(X, Y, np.size(X, 0), theta)
            norm = np.linalg.norm(grad)
        return theta
    
    ''' 以下是测试数据 '''
    
    # 测试用线性函数
    def linear_function(x1, x2, x3):
        result = 1 + 2 * x1 + 3 * x2 + 4 * x3 
        result = result + np.random.rand() # 噪音
        return result
    
    # 计算函数值
    def calculate(X):    
        rowsnumber = np.size(X, axis = 0)    
        Y = [linear_function (X[i, 0], X[i, 1], X[i, 2]) for i in range(0, rowsnumber)]
        return Y
    
    if __name__ == "__main__":
        row_count = 500
        X = np.random.randint(0, 10, (row_count, 3)) # 随机产生row_count个样本
        Y = calculate(X) # 计算标记
    
        X0 = np.ones((row_count, 1)) 
        X = np.hstack((X0, X)) # 补充一列1
    
        theta = gradient_descent(X, Y)
        print('theta is ', theta)
    

    运行结果:theta is [1.41206515 2.00558441 3.0013728 4.00684577]

    上面的迭代方法被称为批量梯度下降法,参考式(4),计算梯度时用到了所有的样本。梯度下降法还有个简化的版本,叫做随机梯度下降法,每次计算梯度时只随机使用一个样本,而不是所有样本,这样可以加快计算速度。将式(4)修改为:

    [frac{partial}{partial heta_n}j( heta_0, heta_1dots heta_n) = ( heta_0x_0^{(k)} + heta_1x_1^{(k)} + dots+ heta_nx_n^{(k)}- y^{(k)})x_n^{(k)} ag{5} ]

    其中:(1 leq k leq m)

    将上面Python代码中的方法gradient替换一下:

    def gradient_sgd(X, Y, m, theta):
        ''' 求theta位置的梯度.
    
        Args:
            X: 样本
            Y: 样本标记
            m: 样本数
            theta: 欲求梯度的位置
        
        Returns:
            gi: theta处函数的梯度值
        '''
        theta_size = np.size(theta)
        g = np.zeros(theta_size)
    
        for i in range(theta_size):   
            random_Index = np.random.randint(1, m + 1)
            gi = ((np.dot(X[random_Index], theta) - Y[random_Index]) * X[random_Index, i])
            g[i] = gi
    
        return g
    

    运行结果:
    theta is [1.43718942 2.00043557 3.00620849 4.00674728]

    感觉像是飞起来。随机梯度下降法肯定没有批量梯度下降法准确,所有还有第三种下降法,叫做小批量梯度下降法,介于批量梯度下降法和随机梯度下降法之间,每次计算梯度使用随机的一小批样本,此处不再code说明。

    正规方程导法

    因为代价函数是个凸函数,那么我们可以对代价函数求导,让其导数等于0的点即为最小值点。

    为方便计算,我们在前面增加了一个值恒等于1的(x_0),这样就把线性函数的偏置项去掉了,参考式(2),重新定义矩阵(X)为:

    [egin{pmatrix} x_0^{(1)} & x_1^{(1)} & x_2^{(1)} & cdots &x_n^{(1)}\ x_0^{(2)} &x_1^{(2)} & x_2^{(2)} & cdots & x_n^{(2)}\ vdots & vdots & vdots & vdots & vdots\ x_0^{(m)} & x_1^{(m)} & x_2^{(m)} & cdots & x_n^{(m)}\ end{pmatrix}]

    代价函数式(3)等价于:

    [J(vec heta)=frac{1}{2m}||Xvec heta - vec{y}||^2 ag{6} ]

    化简式(6):

    [egin{align} J(vec heta)&=frac{1}{2m}||Xvec heta - vec{y}||^2 \ &=frac{1}{2m}(Xvec heta - vec{y})^T(Xvec heta - vec{y}) \ &=frac{1}{2m}(vec heta^TX^T - vec{y}^T)(Xvec heta - vec{y}) \ &=frac{1}{2m}(vec heta^TX^TXvec heta - vec heta^TX^Tvec{y}- vec{y}^TXvec heta + vec{y}^Tvec{y})\ &=frac{1}{2m}(vec heta^TX^TXvec heta - 2vec{y}^TXvec heta + vec{y}^Tvec{y})\ end{align}]

    (vec heta)求导:

    [frac{d}{dvec heta}J(vec heta)=frac{1}{m}(X^TXvec heta-X^Tvec{y}) ]

    令其等于0,得:$$vec heta=(X^TX)^{-1}X^Tvec{y} ag{7}$$

    将上面的Python代码改为:

    # 测试用线性函数
    def linear_function(x1, x2, x3):
        result = 1 + 2 * x1 + 3 * x2 + 4 * x3 
        result = result + np.random.rand() # 噪音
        return result
    
    # 计算函数值
    def calculate(X):    
        rowsnumber = np.size(X, axis = 0)    
        Y = [linear_function (X[i, 0], X[i, 1], X[i, 2]) for i in range(0, rowsnumber)]
        return Y
    
    if __name__ == "__main__":
        row_count = 500
        X = np.random.randint(0, 10, (row_count, 3)) # 随机产生row_count个样本
        Y = calculate(X) # 计算标记
    
        X0 = np.ones((row_count, 1)) 
        X = np.hstack((X0, X)) # 补充一列1
    
        theta = np.dot(np.dot(np.linalg.pinv(np.dot(X.T, X)), X.T), np.array(Y).T)
        print('theta is ', theta)
    

    运行结果:theta is [1.49522638 1.99801209 2.99704438 4.00427252]

    和梯度下降法比较,光速的感觉,那为什么还要用梯度下降法呢?这是因为求矩阵的逆算法复杂度较高,达爷的建议是:如果样本的属性超过一万个,考虑使用梯度下降法。

    调用函数库

    其实我们也可以直接调用类库的,有很多类库可以做回归算法,比如:

    import numpy as np
    from sklearn import linear_model
    
    # 测试用线性函数
    def linear_function(x1, x2, x3):
        result = 1 + 2 * x1 + 3 * x2 + 4 * x3 
        result = result + np.random.rand() # 噪音
        return result
    
    # 计算函数值
    def calculate(X):    
        rowsnumber = np.size(X, axis = 0)    
        Y = [linear_function (X[i, 0], X[i, 1], X[i, 2]) for i in range(0, rowsnumber)]
        return Y
    
    if __name__ == "__main__":
        row_count = 500
        X = np.random.randint(0, 10, (row_count, 3)) # 随机产生row_count个样本
        Y = calculate(X) # 计算标记
    
        regr = linear_model.LinearRegression()
        regr.fit(X, np.array(Y).T)
    
        a, b = regr.coef_, regr.intercept_
        print(a)
        print(b)
    

    运行结果:
    [2.00384674 2.99234723 3.99603084]
    1.5344826581936104

    和我们自己算的差不多吧。还有很多其他的类库可以调用,大叔没有一一去找。可能通常只要调用类库就足够了,不需要我们自己写,不过还是知道原理比较好,遇到问题才好对症下药。

    我是这样理解的:我们能够调用到的常见的(广义)线性回归库,其实内部都是用直接求导法实现的(没有看过源码,猜测是直接求导,如果是梯度下降,不太可能自动算出步长),如果样本的属性比较少,比如少于一万个,调用类库就好,类库肯定比我们大部分人自己写的强,但是当样本属性非常多时,用直接求导法求解速度太慢,这时才需要我们自己写梯度下降代码。

  • 相关阅读:
    geoserver发布地图服务WMTS
    geoserver发布地图服务WMS
    geoserver安装部署步骤
    arcgis api 3.x for js 入门开发系列十四最近设施点路径分析(附源码下载)
    arcgis api 3.x for js 入门开发系列十三地图最短路径分析(附源码下载)
    cesium 之自定义气泡窗口 infoWindow 后续优化篇(附源码下载)
    arcgis api 3.x for js 入门开发系列十二地图打印GP服务(附源码下载)
    arcgis api 3.x for js 入门开发系列十一地图统计图(附源码下载)
    arcgis api 3.x for js 入门开发系列十叠加 SHP 图层(附源码下载)
    arcgis api 3.x for js入门开发系列九热力图效果(附源码下载)
  • 原文地址:https://www.cnblogs.com/zzy0471/p/linear_regression.html
Copyright © 2011-2022 走看看