线性回归算法是机器学习里面最基础、最简单的算法，但是它在实际应用中非常广泛。在接下来的内容，我们将通过一个房价预测问题来引入探讨线性回归算法的数学原理、代价函数、梯度下降算法等。并在分析完原理之后，用python完全实现一遍，并和scikit-learn的机器学习库对比实际的预测效果。

一、基础理论

1、问题描述

        如下图是北京朝阳地区的房价，我们都知道，房价的高低跟很多因素相关，包括：房间面积、楼层位置、房屋位置等等。那要给出一个房间的各种参数，要怎么判断它的房价呢。

地址：https://bj.5i5j.com/ershoufang/chaoyangqu/n1/

        接下来，我们就通过北京朝阳区的房价数据应用线性回归算法进行训练建模，然后对房价进行预测。首先我们需要爬取这些数据，这里偷个懒，使用了网友已经整理好的数据（https://github.com/rieuse/Machine_Learning）。这里爬去的信息包括几房、几厅、房子面积、朝向、楼层位置、装修、单位面积房价和房价。

rooms,halls,size,direction,height,decoration,unit_price,price 2,1,50.26,南,高,简,33824,170 1,1,53.82,南北,低,中,44593,240 2,1,61.01,南北,低,中,54089,330 2,1,52.03,南北,中,简,49971,260 2,1,46.6,南,中,,107296,500 …...

        我们在这里截取rooms、halls、size和price四个字段的数据进行数据训练建模。

        如下图，可以看到：rooms、halls、size分别和price大体是单调递增的关系。

问题：如果通过已知的房价数据进行建模预测后续的房价？

2、模型表示

        房价预测这个问题是一个回归问题，而这里影响房价的因素（变量）有多个，是一个多元回归问题。回归模型正是表示从输入变量到输出变量之间的映射的函数。回归问题分为两个过程，分别是学习和预测，如下图所示。首先是给定一个训练数据集：

        学习系统基于训练数据集构建一个模型（也叫预测函数），即Y=f(X)，对新输入的，预测系统过根据学习的模型Y=f(X)确定相应的输出。

        下面是我们约定的表示方式：

• m代表训练集中样本的数量
• x代表输入变量
• y代表输出变量
• (x,y)代表训练集中的样本
• 代表第i个样本，如上述房价中的一个数据，输入x(1)表示：(rooms：2，halls：1，size：50)，输出y(0)表示：price：170
• 代表第i个样本的第j个分量，举房价数据为参考，表示第一个数据样本中的halls字段，即=1

        房价预测问题，实际上就是通过大量的训练数据来训练模型，最终得到预测函数Y=f(X)。然后我们将要预测的房子的rooms、halls、size作为输入变量，然后通过学习得到的预测函数f(X)得到输出变量为结果。这里我们假设放假的预测函数为：分别代表rooms、halls和size三个变量，所以这里的问题就变成了求解四个参数，找到符合的预测函数。

2、代价函数

        上面，我们假设了预测函数为：。也就是说，我们要求解出四个参数。但是得到参数解以后，我们怎么判定得到的参数是最合适的呢？这里就引用到代价函数的概念。我们知道通过预测函数得到的预测值f(x)和实际的值y是可能存在一定误差的，而我们要做的是就是缩小这些误差。下面我们构建一个代价函数，这个函数是所有建模误差的平方和，即：

        其中。所以，我们的目标就是找出使得代价函数最小的一系列参数变量，代价函数返回的值越小，说明预测函数越拟合数据，当然也可能出现过拟合的情况。

3、梯度下降

3.1 梯度下降算法原理

        梯度下降是一个用来求函数局部最小值的算法，相应地梯度上升算法是用来求函数局部最大值的算法。

        如下图是求通过梯度下降算法求的最小值，可以看到不同的初始值可能导致最终得到不同的局部最小值。

多元线性回归的批量梯度下降算法为：

在上式中α称为学习率，也叫步长，它决定了在梯度下降的迭代过程中，每一步沿梯度负方向前进的长度。以上图的例子，学习率就是这一步所在位置沿着最陡峭的最易下山的位置走的那一步的长度。学习率如果取值太小，会导致收敛变慢，迭代次数变多；取值太大，则会导致波动较大，导致无法收敛。同时，该算法要对所有的θ参数变量进行求导过程。

将上式的代价函数展开，即：

求导后得到：

这里，对求导后，分别为：

下面，我们对梯度下降算法做个整理：

3.2  梯度下降算法的矩阵方式描述

在3.1中，我们推导得到了梯度下降算法的代数描述方式。

下面我们用更简洁的矩阵方式来描述梯度下降算法。

首先约定如下几个矩阵形式：

令：

综合起来：

所以，采用矩阵形式表示梯度下降后，θ更新步骤如下，通过这个公式，我们就可以通过pandas和numpy更简洁地实现代码。

备注：这里θ是nx1矩阵，x是mxn矩阵，y是mx1矩阵

3.4、代码实现

        这一节，我们来用python手动实现线性回归算法，通过上面几节的分析，我们可以总结出几个关键的点

（1）预测函数

        备注：上述公式x是nx1矩阵，θ是nx1矩阵，即代表一个数据实例的计算结果

        上式是对一个数据实例的计算结果，如果要对所有m个训练集实例进行矩阵运算，可以采用下述公式：

        备注：上述公式x是mxn矩阵，θ是nx1矩阵，即对所有训练集的数据实例的计算结果

        所以实现代码如下：

def model(theta, X):
    """
    预测函数
    :param theta: 参数矩阵(nx1)
    :param X:输入变量(矩阵mxn)
    :return: 输出变量, mx1
    """
    return X.dot(theta)

（2）代价函数

def cost(m, theta, X, y):
    """
    代价函数
    :param m: 训练集个数
    :param theta: 参数矩阵(mx1)
    :param X: 输入变量矩阵(mx1)
    :param y: 输出变量矩阵(mx1)
    :return:
    """
    ele = model(theta, X) - y
    item = ele ** 2
    item_sum = np.sum(item)
    return item_sum / (2 * m)

（3）梯度下降和θ参数更新

在前几节中，我们推到了上述公式，如下：

def derivative(m, theta, X, y, n):
    """
    对各个theta参数求导
    :param m: 训练集个数
    :param theta: 参数矩阵(nx1)
    :param X: 输入变量矩阵(mxn)
    :param y: 输出变量矩阵(mx1)
    :param n: 𝜃变量个数𝜃0, 𝜃1, 𝜃2, ..., 𝜃n-1
    :return: nx1 矩阵
    """
    grad_theta = []
    for j in range(n):
        # 对所有𝜃j变量求导
        grad = (model(theta, X) - y).T.dot(X[:, j])
        grad_sum = np.sum(grad)
        grad_theta.append(grad_sum / m)
    return np.array(grad_theta).reshape(n, 1)

 

当然如果采用向量形式，代码会更简洁，在上面的章节中，我们推导出了如下的向量表达形式：

def derivative(m, theta, X, y, n):
    """
    对各个theta参数求导（向量形式）
    :param m: 训练集个数
    :param theta: 参数矩阵(nx1)
    :param X: 输入变量矩阵(mxn)
    :param y: 输出变量矩阵(mx1)
    :param n: 𝜃变量个数𝜃0, 𝜃1, 𝜃2, ..., 𝜃n-1
    :return: nx1 矩阵
    """
    grad_sum = X.T.dot(model(theta, X) - y)
    return grad_sum / m

对比下两个代码，明显向量形式的更简洁、更容易理解。

def gradient(grade_theta, theta, alpha):
    """
    𝜃更新
    :param grade_theta:求导后的𝜃矩阵(nx1)
    :param theta: 更新后的𝜃矩阵(nx1)
    :param alpha: 学习率
    :return: 返回更新后的𝜃矩阵(nx1)
    """
    return theta - alpha * grade_theta

（4）停止迭代策略（代价函数小于阈值或设定迭代次数）

def stop_strage(cost, cost_update, threshold):
    """
    停止迭代条件
    :param cost:
    :param cost_update:
    :param threshold:
    :return:
    """
    return (cost >= cost_update) and (cost - cost_update < threshold)

（5）线性回归

def linear_regression(X, y):
    """
    线性回归模型
    :param X: 输入变量矩阵(mxn)
    :param y: 输出变量矩阵(mx1)
    :param threshold:
    :return:
    """
    start = time.clock()
    m = X.shape[0] # 样本个数
    n = X.shape[1]
    theta = np.zeros(n).reshape(n, 1) # 设置权重参数的初始值
    y = y.reshape(m, 1)
    print "theta:", theta.shape, ", X: ", X.shape, "y:", y.shape

    cost_record = [] # 记录代价函数的值
    alpha = 0.0001 # 学习率
    threshold = 0.01
    cost_val = cost(m, theta, X, y)

    cost_record.append(cost_val)

    iters = 0
    while True:
        grad = derivative(m, theta, X, y, n)
        # 权重参数更新
        theta = gradient(grad, theta, alpha)
        cost_update = cost(m, theta, X, y)
        if stop_strage(cost_val, cost_update, threshold):
            break
        cost_val = cost_update
        cost_record.append(cost_val)
        iters += 1
    end = time.clock()
    print("cost time: %f s" % (end - start))

    return cost_record, theta, iters

（6）使用训练线性回归模型预测房价

df = pd.read_csv('house_price.csv')

print df['rooms'].unique()
print df['halls'].unique()

# 去除脏数据
keep_index = (df['rooms'] != '多') & (df['halls'] != '多')
df = df[keep_index]

# 转换数据类型
df['rooms'] = df['rooms'].astype('int')
df['halls'] = df['halls'].astype('int')

# 截取rooms、halls、size和price
X = df[['rooms', 'halls', 'size']]
y = df['price']

cost_record, theta, iters = linear_regression(X.values, y.values)
print iters
print theta
print cost_record

# 预测数据
# rooms,halls,size,price
# 1,1,49,326
# 2,1,56.45,315
# 2,1,58.3,466
# 3,1,92.33,1061
# 2,1,76.88,408

print model(theta, np.array([1, 1, 49]))
print model(theta, np.array([2, 1, 56.45]))
print model(theta, np.array([2, 1, 58.3]))
print model(theta, np.array([3, 1, 92.33]))
print model(theta, np.array([2, 1, 76.88]))

预测得到的房价如下：最终一共迭代iters=11652次，得到的theta=(13.13249986, 12.4409009, 6.59238776)。

[ 348.60040077]
[ 410.84618941]
[ 423.04210675]
[ 660.51356193]
[ 545.52867125]

可以看到部分数据还是预测的比较准的。

下面，我们再来看下代价函数返回的值的收敛趋势图：

# 第一个数据太大，剔除，方便观察对比
plt.plot(np.arange(len(cost_record)-1), cost_record[1:], 'r')
plt.show()

二、进阶

1、梯度下降算法优化

（1）特征缩放（标准化）

在我们上面使用的房价数据里使用了多维特征，如：rooms、halls和size。但是这几个特征的取值范围各不一样。如下：

这种多维特征取值差异大的会影响梯度下降的收敛效果。例如下面吴恩达讲义里介绍的房价的例子，房间数和房子面积两个特征取值相差很大，成左图的一个椭圆形，这会导致梯度下降是来回振荡，影响收敛速度。

解决的办法是，将所有的特征都归一化，让所有特征的取值都缩放到-1到1之间，如下图的右侧所示。可见，这是梯度下降就很稳定，不会来回振荡，这也使得梯度下降算法收敛速度加快。

常见的特征缩放（标准化）方法有

z-score：

min-max：

三、scikit-learn实践

在第一部分，我们通过分析线性回归算法原理，自己动手实现了一个线性回归算法模型来预测房价。这一节内容，我们通过scikit-learn第三方算法库的线性回归工具来预测房价，看看跟我们预测的房价是否相似。

备注：具体的函数使用方法参见scikit-learn的文档：http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html

同样采用的训练数据是一样的。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

df = pd.read_csv('house_price.csv')

# 去除脏数据
keep_index = (df['rooms'] != '多') & (df['halls'] != '多')
df = df[keep_index]

# 转换数据类型
df['rooms'] = df['rooms'].astype('int')
df['halls'] = df['halls'].astype('int')

# 截取rooms、halls、size和price
X = df[['rooms', 'halls', 'size']]
y = df['price']

print "\n======================scikit-learn第三方库的线性回归算法==========================="
#测试集和训练集
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

# 训练
linreg = LinearRegression()
linreg.fit(X_train, y_train)

#结果
print "theta参数："
print linreg.coef_

# 预测房价
print "输入：1,1,49 实际值：326,预测值：", linreg.predict(np.array([1, 1, 49]).reshape(1,3))
print "输入：2,1,56.45 实际值：315,预测值：", linreg.predict(np.array([2, 1, 56.45]).reshape(1,3))
print "输入：2,1,58.3 实际值：466,预测值：", linreg.predict(np.array([2, 1, 58.3]).reshape(1,3))
print "输入：3,1,92.33 实际值：1061,预测值：", linreg.predict(np.array([3, 1, 92.33]).reshape(1,3))
print "输入：2,1,76.88 实际值：408,预测值：", linreg.predict(np.array([2, 1, 76.88]).reshape(1,3))

输出：

theta参数： [ 16.94375062  61.3768412 6.59171506] 输入：1,1,49 实际值：326,预测值： [ 329.56288307] 输入：2,1,56.45 实际值：315,预测值： [ 395.61491088] 输入：2,1,58.3 实际值：466,预测值： [ 407.80958374] 输入：3,1,92.33 实际值：1061,预测值： [ 649.06939779] 输入：2,1,76.88 实际值：408,预测值： [ 530.28364952]

 

跟我们自己实现的线性回归算法预测出来的值还是很相似的，这说明我们的算法还是有效的。

	rooms	halls	size	price	自己实现的线性回归算法	采用scikit-learn的线性回归算法
1	1	1	49	326	348.60040077	329.56288307
2	2	1	56.45	315	410.84618941	395.61491088
3	2	1	58.3	466	423.04210675	407.80958374
4	3	1	92.33	1061	660.51356193	649.06939779
5	2	1	76.88	408	545.52867125	530.28364952