机器学习之线性回归：算法兑现为python代码

前面三天推送机器学习线性回归算法之最小二乘法，从假设到原理，详细分析了直接求解和梯度下降两种算法，接下来手动编写python代码实现线性回归的算法吧。

01 数据预处理

在拿到一个数据集后，往往需要经过漫长的预处理过程，不要忽视这个看似与建立模型，求解模型无关的步骤，它其实非常重要的，为后续工作做好准备的一步。现在这节的重点不是在论述预处理的方法，所以在此不详细介绍预处理的过程，而是重点例子模拟线性回归最小二乘法的两个求解方法。

获得了数据集后，经过预处理后得到的数据前10条展示如下，其中第一列为房屋的面积，第二列为房屋使用年限，第三列为房屋的价值，是标签值，这些值都已经经过预处理了。

房屋面积使用年限价值

[[ 0.35291809, 0.16468428, 0.35774628],

[-0.55106013, -0.10981663, 0.25468008],

[-0.65439632, -0.71406955, 0.1061582 ],

[-0.19790689, 0.61536205, 0.43122894],

[-0.00171825, 0.66827656, 0.44198075],

[-0.2739687 , -1.16342739, 0.01195186],

[ 0.11592071, -0.18320789, 0.29397728],

[-0.02707248, -0.53269863, 0.21784183],

[ 0.7321352 , 0.27868019, 0.42643361],

[-0.76680149, -0.89838545, 0.06411818]]

下面用直接求法和梯度下降法求解线性回归。

首先介绍下使用的库：

'导入numpy库'

import numpy as np

'导入pyplot'

import matplotlib.pyplot as plt

'导入时间模块'

import time

本次模拟取数据集的前100条数据进行迭代计算，即 m = 100 。

做一个偏移量和2个特征的组合，这样与前面推送的理论部分衔接在一起，组合的代码

如下所示：

'偏移量 b shape=(100,1)'

b = np.array([1])

b=b.repeat(100)

'将偏移量与2个特征值组合 shape = (100,3)'

X = np.column_stack((b,X))

02 直接求解参数

我们知道当我们建立线性回归的模型时，因为是线性的，并且误差项满足高斯分布，此时借助最大似然估计可以直接拿到权重参数的计算公式，如果想看下理论部分，请参考直接求解：

'这是一个求矩阵的逆所用到的模块'

from numpy.linalg import linalg as la

xtx = X.transpose().dot(X)

xtx = la.inv(xtx)

'直接求出参数'

theta = xtx.dot(X.transpose()).dot(y)

这个解法很简单，直接套用公式，求出权重参数如下：

array([0.29348374, 0.10224818, 0.19596799])

即偏移量为 0.29，第一个特征的权重参数为0.10，第二个特征的权重参数为0.195 。

下面用梯度下降法求解，这才是我们论述的重点，这个思路与机器学习的其他算法，比如逻辑回归等思路是一致的，因此有必要好好研究下。

03 梯度下降求参数

梯度下降的详细介绍，请参考梯度下降求解权重参数部分，下面我们论述如何由理论兑现为代码。

首先列举梯度下降的思路步骤，采取线性回归模型，求出代价函数，进而求出梯度，求偏导是重要的一步，然后设定一个学习率迭代参数，当与前一时步的代价函数与当前的代价函数的差小于阈值时，计算结束，如下思路：

'model' 建立的线性回归模型
'cost' 代价函数
'gradient' 梯度公式
'theta update' 参数更新公式
'stop stratege' 迭代停止策略：代价函数小于阈值法

下面分别写出以上五步的具体实现代码，

'model'
def model(theta,X):
    theta = np.array(theta)
    return X.dot(theta)
'cost'
def cost(m,theta,X,y):
    'print(theta)'
    ele = y - model(theta,X)
    item = ele**2
    item_sum = np.sum(item)
    return item_sum/2/m
'gradient'
def gradient(m,theta,X,y,cols):
    grad_theta = []
    for j in range(cols):
        grad = (y-model(theta,X)).dot(X[:,j])
        grad_sum = np.sum(grad)    
        grad_theta.append(-grad_sum/m)
    return np.array(grad_theta)
'theta update'
def theta_update(grad_theta,theta,sigma):
    return theta - sigma * grad_theta
'stop stratege'
def stop_stratege(cost,cost_update,threshold):
    return cost-cost_update < threshold
'OLS algorithm'
def OLS(X,y,threshold):
    start = time.clock()
    '样本个数'
    m=100
    '设置权重参数的初始值'
    theta = [0,0,0]
    '迭代步数'
    iters = 0;
    '记录代价函数的值'
    cost_record=[]
    '学习率'
    sigma = 0.0001
    cost_val = cost(m,theta,X,y)
    cost_record.append(cost_val)
    while True:
        grad = gradient(m,theta,X,y,3)
        '参数更新'
        theta = theta_update(grad,theta,sigma)
        cost_update = cost(m,theta,X,y)
        if stop_stratege(cost_val,cost_update,threshold):
            break
        iters=iters+1
        cost_val = cost_update
        cost_record.append(cost_val)
    end = time.clock()
    print("OLS convergence duration: %f s" % (end - start))
    return cost_record, iters,theta

结果显示经过，OLS梯度下降经过如下时间得到初步收敛，OLS convergence duration: 7.456927 s，经过3万多个时步迭代，每个时步计算代价函数的取值，如下图所示：

收敛时，得到的权重参数为：

array([ 0.29921652, 0.09754371, 0.1867609 ])

可以看到梯度下降得到的权重参数与直接求出法得出的基本相似，这其中的误差是因为没有进一步再迭代。

04 总结

以上就是最小二乘法的两种解法的代码实现，至此我们已经将线性回归算法的最基本的OLS从理论，假设，到现在的代码兑现都阐述完了。让我们看一下远边的大海，和巍峨的高山，放松一下吧！

然而，有些数据集的某两列或多列存在强相关性，当面对这样的数据集，OLS还能胜任吗？如果不能胜任，这其中的原因又是什么呢？

请看明天的推送，OLS算法的缺陷及原理。