[机器学习] 房屋价格预测(补充)

前言

在之前有许多关于线性回归-房屋价格预测的优秀文章, 这篇blog更适合看了吴恩达机器学习，想通过Python实现线性回归。我试着写了代码，在过程中遇到一些问题，有的是公式问题，有的是关于numpy库的问题，有的blog在这些方面不太完善，我在这里会一一列举出来。

项目地址

暂无

参考连接

吴恩达机器学习 线性回归

前导知识

需要知道关于假设函数
$$H(x)={\Theta }_0+{\Theta }_{1}x$$
需要知道关于代价函数
$$J({\Theta }_0,{\Theta }_1)=\frac{{\sum }_{i=1}^m[h(x^{(i)})-y^{(i)}{]}^2}{2m}$$
需要知道关于多元代价函数的偏导数
$$J({\Theta }_j{)}^{\prime }=\frac{{\sum }_{i=1}^m(h(x^{(i)})-y^{(i)})x^{(i)}}{m}$$

需要知道关于均值归一化
$$X^{(i)}=\frac{X^{(i)}-\mu }{X_{max}-X_{min}}$$

还需要知道一些numpy的基本用法, 以及matplotlib的简单使用。

正文

一、加载数据集

def load_file(filepath):
    """
        加载文件, 数据类型默认float64
    :param filepath: 文件
    :return: m * n 的 array, m 为数据条目, n 为特征条目
    """
    return np.loadtxt(filepath, delimiter=',', dtype=np.float64)

load_file函数返回np.array数据类型的数据, 例如:
$$\left[\begin{array}{cc}2104& 399900\\ 1600& 329900\\ 1416& 232000\end{array}\right]$$
以 np.array 表示为:
$$\begin{gathered} [[2104,399900] \\ [1600,329900] \\ [1416,232000]] \end{gathered}$$
是一个$3\ast 2$的矩阵

二、绘制数据

import numpy as np
from matplotlib import pyplot as plt

def plot_house_price(data, theta = None):
    """
        绘制房屋面积与价格的离散图
    :param data: 为 np.array 类型 含两列数据, 第0列指面积, 第1列指房屋价格
    :param theta: [Θ₀, Θ₁] 拟合函数系数
    """
    # 设置 x and y 轴的标签 与 该图形的标题
    plt.xlabel("area")
    plt.ylabel("price")
    plt.title("Houses sale instance")
    # 绘制圆点, data[i][0] 为面积, data[i][1] 为价格
    plt.scatter(data[:, 0], data[:, 1], s=20)

    # 如果传入theta列表, 表示需要绘制拟合函数
    if theta is not None:
        # 设置显示范围
        x_axis = np.linspace(-np.max(data[:, 0]) * 1.2, np.max(data[:, 0]) * 1.2)
        # h(x) = Θ₀ + Θ₁ x
        y_axis = theta[0] + theta[1] * x_axis
        # 绘制拟合函数
        plt.plot(x_axis, y_axis, "g-")

    # 显示
    plt.show()

其中, 一些绘图的操作便不再赘述。值得注意的是$data[:,0]$，这是一个切片操作，逗号前后分别是行索引和列索引的切片，那么就这句意为“取矩阵的第0列(所有行)”。

三、对数据集归一化

归一化的作用是为了使函数快速收敛。如果对于差异度很大的数据集，如果不归一化，使用梯度下降算法训练模型时，可能导致函数不收敛，也就无法训练出模型。（本例即是）

def mean_normalization(data):
    """
        data 为 np.array格式 含两列数据, 第0列指面积, 第1列指房屋价格
        均值归一化
    :return:
    """

    # col_mean 即 列的均值; col_std 即 列的标准差;
    # 参数 0 代表列; 若求行均值或标准差 传入参数 1 即可
    col_mean = np.mean(data, 0)
    col_max = np.max(data, 0)
    col_min = np.min(data, 0)

    # shape 返回 维度的维数 即
    # l = [[1, 2, 3], [4, 5, 6]] => shape[0] = 2, shape[1] = 3
    for i in range(data.shape[1]):
        """
            shape[1] 返回列数
             Xᵢ = (Xᵢ - μ) / (Xₘₐₓ - Xₘᵢₙ)
            其中, Xᵢ 为第i个特征, μ 是当前特征的均值, (Xₘₐₓ - Xₘᵢₙ) 是当前特征的范围
        """

        # 让 当前特征中每个特征值都执行归一化操作
        # 均值归一化, 一个特征即矩阵中一列
        data[:, i] = (data[:, i] - col_mean[i]) / (col_max[i] - col_min[i])

# 这段代码，虽无必要说明，但还是贴到这里
# 加载资源，然后绘制房屋价格点图，然后执行归一化后，再绘制
if __name__ == '__main__':
    data = load_file("data.txt")
    plot_house_price(data)
    print(data)
    mean_normalization(data)
    plot_house_price(data)
    print(data)

对比没有归一化之前的图像, 两者差别不大。

番外：绘制J函数

停停… 在对数据集归一化后，就可以执行梯度下降算法，解出Θ值。但是，我在做的过程中，对代价函数J非常感兴趣，所以耗费一点时间绘制出来。这里附带一个推导：

先说结论，
$$J({\Theta }_0,{\Theta }_1)=\frac{{\sum }_{i=1}^m[h(x^i)-y^i{]}^2}{2m}=\frac{K(x\Theta -y{)}^2}{2m}=\frac{K{A}^2}{2m}$$
推导如下，
设$x$为列向量，其值有$$\left[\begin{array}{c}{x}^{(1)}\\ x^{(2)}\\ x^{(3)}\\ x^{(4)}\\ \dots \\ x^{(m)}\end{array}\right]$$
设$y$同为列向量，其值有$$\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ y^{(3)}\\ y^{(4)}\\ \dots \\ y^{(m)}\end{array}\right]$$
系数矩阵$\Theta$为，
$$\left[\begin{array}{c}{\Theta }_0\\ {\Theta }_1\end{array}\right]$$
我们知道假设函数$h(x)={\Theta }_0+{\Theta }_{1}x$，那么$h(x^i)={\Theta }_0+{\Theta }_1{x}^i$，有$h(x^i)$的矩阵形式：
$$\left[\begin{array}{c}{\Theta }_0+{\Theta }_1{x}^{(1)}\\ {\Theta }_0+{\Theta }_1{x}^{(2)}\\ {\Theta }_0+{\Theta }_1{x}^{(3)}\\ {\Theta }_0+{\Theta }_1{x}^{(4)}\\ \dots \\ {\Theta }_0+{\Theta }_1{x}^{(m)}\end{array}\right]$$
所以，为了运算方便，只需给$x$加一列$1$，使得$x^{\prime }$ 为
$$\left[\begin{array}{cc}1& x^{(1)}\\ 1& x^{(2)}\\ 1& x^{(3)}\\ 1& x^{(4)}\\ \dots & \dots \\ 1& x^{(m)}\end{array}\right]$$，那么，可以推出$h(x^i)=x^{\prime }\Theta$(注意,$\Theta$是系数矩阵)，则有矩阵$A=h(x^i)-y^i=x^{\prime }\Theta -y$，得出矩阵$B$，$B=A^2=[h(x^i)-y^i{]}^2$。
如何求解${\sum }_{i=1}^m$呢？这里, 可以利用矩阵乘法的性质：
已知矩阵B为
$$\left[\begin{array}{c}{A}^{2(1)}\\ A^{2(2)}\\ A^{2(3)}\\ A^{2(4)}\\ \dots \\ A^{2(m)}\end{array}\right]$$
,取矩阵$K$为
$$\left[\begin{array}{cccc}1& 1& \dots & 1\end{array}\right]$$
, 大小为$1\ast m$, 矩阵$B$的大小为$m\ast 1$，则$K\ast B=S$, $S$即$$\sum _{i=1}^m[h(x^i)-y^i{]}^2$$

def plot_J_function(sr_x, sr_y):
    """
        绘制代价函数J
    :param sr_x: 房屋面积, 归一化后的数据
    :param sr_y: 房屋价格, 同归一化后的数据
    :return: 
    """
    
    # 样本数量
    m = sr_x.shape[0]

    # 扩展sr_x
    sr_x = np.column_stack((np.ones(len(sr_x)), sr_x))
    # 全 1 行向量
    K = np.array([np.ones(sr_x.shape[0])])

    # X and Y轴 刻度
    delta = 1
    # X, Y 的显示范围
    X = np.arange(-10, 10, delta)
    Y = np.arange(-10, 10, delta)

    Z = []
    for t0 in X:
        Zl = []
        for t1 in Y:
            theta = np.array([[t0, t1]]).T
            J = np.dot(K,
                       np.power(np.dot(sr_x, theta) - sr_y, 2)
                ) / (2 * m)
            Zl.append(J[0][0])
        Z.append(Zl)
    Z = np.array(Z)

    # 扩充, 此时X,Y分别是一个矩阵, 内含所有可能的取样点
    X, Y = np.meshgrid(X, Y)
    # 显示
    ax = plt.axes(projection='3d')
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.cm.coolwarm, edgecolor='none')
    ax.set_xlabel('theta_0')
    ax.set_ylabel('theta_1')
    ax.set_zlabel('J')
    plt.show()

if __name__ == '__main__':
    data = load_file("data.txt")
    plot_house_price(data)
    mean_normalization(data)
    plot_house_price(data)
    x = np.array([data[:, 0]]).T
    y = np.array([data[:, 1]]).T
    plot_J_function(x, y)

四、梯度下降求解Θ

要找到$minmizeJ({\Theta }_0,{\Theta }_1)$的极小值点，只需要求出$J$的各项偏导, 再令其偏导数 = 0，就可以解出${\Theta }_0,{\Theta }_1,\dots ,{\Theta }_n$。

有J的偏导数，
$$J({\Theta }_j{)}^{\prime }=\frac{{\sum }_{i=1}^m(h(x^{(i)})-y^{(i)})x^{(i)}}{m}$$，规定 $x^{(0)}=1$, 则有
$$\begin{gathered} J({\Theta }_0{)}^{\prime }=\frac{{\sum }_{i=1}^m(h(x^{(i)})-y^{(i)})}{m} \\ J({\Theta }_1{)}^{\prime }=\frac{{\sum }_{i=1}^m(h(x^{(i)})-y^{(i)})x^{(i)}}{m} \end{gathered}$$在绘制J函数部分有讲到，矩阵$A=h(x^i)-y^i=x^{\prime }\Theta -y$，$K\ast B=S$，而$K$是全$1$的行向量，其实，在这里，只需将$K$替换成$x^T$即可，所以有
$$\sum _{i=1}^m(h(x^{(i)})-y^{(i)})x^{(i)}=x^T(x^{\prime }\Theta -y)$$，这部分可以自己尝试着证明一下。
所以，
$${\Theta }_j={\Theta }_j-\frac{\alpha }{m}J({\Theta }_j{)}^{\prime }$$ 随着迭代次数的增多, $J({\Theta }_j{)}^{\prime }慢慢趋向0$，最后${\Theta }_j$便稳定下来。

def gradient_descent(sr_x, sr_y, theta, alpha, iter_steps):
    """
        梯度下降
    :param sr_x: 列向量
    :param sr_y: 列向量
    :param theta: 列向量
    :param alpha: 学习率
    :param iter_steps: 迭代次数
    :return: theta
    """
    theta0 = theta[0][0]
    theta1 = theta[1][0]

    # 样本数量
    m = sr_x.shape[0]
    # 扩展sr_x
    X = np.column_stack((np.ones(len(sr_x)), sr_x))
    xT = sr_x.T
    K = np.array([np.ones(sr_x.shape[0])])

    for i in range(iter_steps):
        # XΘ-Y
        hypothesis = np.dot(X, theta)
        loss = hypothesis - sr_y

        # update theta0
        temp0 = theta[0][0] - (alpha / m) * np.dot(K, loss)
        # update theta1
        temp1 = theta[1][0] - (alpha / m) * np.dot(xT, loss)
        theta[0][0] = temp0
        theta[1][0] = temp1

    return theta

if __name__ == '__main__':
    data = load_file("data.txt")
    plot_house_price(data)
    mean_normalization(data)
    plot_house_price(data)
    theta = np.array([[1., 1.]]).T
    theta = gradient_descent(x, y, theta, 0.05, 500)
    print(theta)
    plot_house_price(data, [theta[0][0], theta[1][0]])

小结

本来还有一些东西要讲，但是我换了个思路写代码，那些小疑问也就不需要再说了。例如在梯度下降时，别人的代码是这样的：

    xTrains = x.transpose()
    for i in range(0, maxIterations):
        hypothesis = np.dot(x, theta)
        loss = hypothesis - y
        gradient = np.dot(xTrains, loss) / m
        theta = theta - alpha * gradient
    return theta

与我的差别不大，利用的是一维数组点乘（不知道起什么名…），在我第一次看到这段代码，并运行出结果时感觉非常惊讶，可能是不太熟悉numpy的原因，这段代码中的theta是[1, 1] – 一个list。操作有点迷，但是理解之后便觉得非常好用。
解释 - numpy数组与矩阵乘法
我的代码则按照吴恩达老师的视频讲解那样，先更新theta0，再更新theta1，最后得出结果。

第一次写线性回归时，还有个疑问（没认真听的后果），theta的初始值该为什么？，这个问题想明白也会觉得‘不过如此’。
当然还有一个关于np.array的问题也提点一句，np.array([1, 2])返回的是一维数组，而不是向量；np.array([[1, 2]])这才是行向量。可以用shape属性输出就明白了。

结疑

无

标签

***机器学习，吴恩达视频，房屋价格预测 ***