基于线性回归（Linear Regression）的房屋价格预测

线性回归是统计学中的一种基本预测模型，用于估计因变量（响应变量）和一个或多个自变量（解释变量）之间的关系。线性回归模型假设这些变量之间存在线性关系。根据自变量的数量，线性回归可以分为简单线性回归（一个自变量）和多元线性回归（多个自变量）。

梯度下降是一种优化算法，用于最小化一个函数，通常用于机器学习中的参数优化问题。其核心思想是：通过迭代地调整参数，沿着目标函数（通常是损失函数）梯度下降的方向逐步逼近最小值。

线性回归的损失函数是凸的，可以使用梯度下降法来拟合线性回归模型的参数，涉及以下几个步骤：

1. 初始化参数。随机初始化线性回归模型的参数（权重和偏置）。对于简单线性回归，通常是初始化权重$\theta$和偏置$b$（有时也将偏置包含在权重中）。

2. 定义损失函数。通常使用均方误差（Mean Squared Error, MSE）作为损失函数。对于给定的训练数据集$(x_i,y_i)$，损失函数定义为：$J(\theta )=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x_i)-y_i{)}^2$。其中$h_{\theta }(x_i)={\theta }^T{x}_i+b$是模型的预测值，$m$是样本数量。

3. 计算梯度。计算损失函数关于每个参数的梯度。对于每个参数${\theta }_j$，梯度计算为：$\frac{\partial J(\theta )}{\partial {\theta }_j}=\frac{2}{m}\sum _{i=1}^m(h_{\theta }(x_i)-y_i)x_{ij}$；对于偏置项$b$，梯度为：$\frac{\partial J(\theta )}{\partial b}=\frac{2}{m}\sum _{i=1}^m(h_{\theta }(x_i)-y_i)$

4. 参数更新：使用梯度下降法更新参数。对于每个参数${\theta }_j$和偏置$b$，更新公式为：${\theta }_j:={\theta }_j-\alpha \frac{\partial J(\theta )}{\partial {\theta }_j}$；$b:=b-\alpha \frac{\partial J(\theta )}{\partial b}$
   其中$\alpha$是学习率，决定了每次更新的步长。

5. 重复迭代。重复步骤3和步骤4，直到损失函数收敛到一个较小的值，或者达到预设的迭代次数。

6. 终止条件。可以根据损失函数的变化情况（如变化小于某个阈值）或达到最大迭代次数来终止训练。

7. 模型评估：使用测试数据集评估模型的性能，计算预测值与真实值之间的误差。

一元线性回归

数据集介绍：吴恩达作业（ex1data1.txt）

6.1101,17.592
5.5277,9.1302
8.5186,13.662
7.0032,11.854
5.8598,6.8233
8.3829,11.886
7.4764,4.3483
8.5781,12
6.4862,6.5987
5.0546,3.8166
5.7107,3.2522
14.164,15.505
5.734,3.1551
8.4084,7.2258
5.6407,0.71618
5.3794,3.5129
6.3654,5.3048
5.1301,0.56077
6.4296,3.6518
7.0708,5.3893
6.1891,3.1386
20.27,21.767
5.4901,4.263
6.3261,5.1875
5.5649,3.0825
18.945,22.638
12.828,13.501
10.957,7.0467
13.176,14.692
22.203,24.147
5.2524,-1.22
6.5894,5.9966
9.2482,12.134
5.8918,1.8495
8.2111,6.5426
7.9334,4.5623
8.0959,4.1164
5.6063,3.3928
12.836,10.117
6.3534,5.4974
5.4069,0.55657
6.8825,3.9115
11.708,5.3854
5.7737,2.4406
7.8247,6.7318
7.0931,1.0463
5.0702,5.1337
5.8014,1.844
11.7,8.0043
5.5416,1.0179
7.5402,6.7504
5.3077,1.8396
7.4239,4.2885
7.6031,4.9981
6.3328,1.4233
6.3589,-1.4211
6.2742,2.4756
5.6397,4.6042
9.3102,3.9624
9.4536,5.4141
8.8254,5.1694
5.1793,-0.74279
21.279,17.929
14.908,12.054
18.959,17.054
7.2182,4.8852
8.2951,5.7442
10.236,7.7754
5.4994,1.0173
20.341,20.992
10.136,6.6799
7.3345,4.0259
6.0062,1.2784
7.2259,3.3411
5.0269,-2.6807
6.5479,0.29678
7.5386,3.8845
5.0365,5.7014
10.274,6.7526
5.1077,2.0576
5.7292,0.47953
5.1884,0.20421
6.3557,0.67861
9.7687,7.5435
6.5159,5.3436
8.5172,4.2415
9.1802,6.7981
6.002,0.92695
5.5204,0.152
5.0594,2.8214
5.7077,1.8451
7.6366,4.2959
5.8707,7.2029
5.3054,1.9869
8.2934,0.14454
13.394,9.0551
5.4369,0.61705

餐馆老板想在不同城市开分店，并收集了这些城市的 人口数和当地餐馆的利润作为训练样本保存在 ex1data1.txt。该数据集为97行×2列，第一列为人口数（万），第二列为餐馆利润（万元）。请利用线性回归，根据城市的人口数据预测其利润。

加载数据

导入库：

import numpy as np  			# 科学计算库，处理多维数组，进行数据分析
import pandas as pd				# 基于 NumPy 的一种解决数据分析任务工具
import matplotlib.pyplot as plt # 提供一个类似 Matlab 的绘图框架

导入数据集：

data = np.loadtxt("./resources/ex1data1.txt", delimiter=',', dtype=np.float64)

数据级 D = { ( x ( i ) , y ( i ) ) } i = 1 N = 97 \mathcal{D}=\{{(x^{(i)}, y^{(i)})} \}_{i=1}^{N=97} D={(x(i),y(i))}i=1N=97​，$x^{(i)}$为第$i$个样本，即第$i$行的人口数，$N$为样本数量。

假设函数

已知单变量线性回归的假设函数为$f(x)=wx+b$，进一步向量乘法的形式，表示起来更简单，易于编程计算，这里构建增广权重向量$w=\left[\begin{array}{c}w\\ b\end{array}\right]$，增广特征向量$x=\left[\begin{array}{c}x\\ 1\end{array}\right]$：
KaTeX parse error: {align} can be used only in display mode.

获取数据集中的输入向量$X$，需要在数据集加入$x_0=1$的一列：

data = np.loadtxt("ex1data1.txt", delimiter=",", dtype=np.float64)
y = data[:, 1:2]
X = np.insert(data[:, 0:1], obj=1, values=1, axis=1)

初始化参数$w=\left[\begin{array}{c}0\\ 0\end{array}\right]$为2行1列

w = np.zeros(X.shape[1]).reshape(2, 1)

损失函数（代价函数）

令$X=\left[\begin{array}{cc}{x}^{(1)}& 1\\ x^{(2)}& 1\\ ⋮& ⋮\\ x^{(N)}& 1\end{array}\right]$，$y=\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(N)}\end{array}\right]$，则损失函数$\mathcal{L}(w)$可以表示为：

KaTeX parse error: {align} can be used only in display mode.

def loss_func(w, X, y):
    return np.sum(np.square(X.dot(w) - y)) / 2 * X.shape[0]

梯度下降函数

首先对$\mathcal{L}(w)$求导，KaTeX parse error: {align} can be used only in display mode.，这里用到了向量函数及其导数公式：$\frac{\partial Ax}{\partial x}=A^{\mathrm{T}}$。故梯度下降更新公式：
KaTeX parse error: {align} can be used only in display mode.

def grad(w, X, y):
    return (X.T @ (X.dot(w) - y)) / X.shape[0]

def grad_descent(w, X, y, alpha, maxIter):
    loss = []
    for i in range(maxIter):
        w = w - alpha * grad(w, X, y)
        cost = loss_func(w, X, y)
        loss.append(cost)
    return w, loss

完整代码

import numpy as np
from matplotlib import pyplot as plt


def loss_func(w, X, y):
    return np.sum(np.square(X.dot(w) - y)) / 2 * X.shape[0]


def grad(w, X, y):
    return (X.T @ (X.dot(w) - y)) / X.shape[0]


def grad_descent(w, X, y, alpha, maxIter):
    loss = []
    for i in range(maxIter):
        w = w - alpha * grad(w, X, y)
        cost = loss_func(w, X, y)
        loss.append(cost)
    return w, loss


def paint(loss, data):
    fig = plt.figure(figsize=(10, 4))
    # 画第1个图：折线图
    ax1 = fig.add_subplot(121)
    ax1.plot([i for i in range(maxIter)], loss, color="#FF0000")
    ax1.set_title("loss function")
    # 画第2个图：散点图
    ax2 = fig.add_subplot(122)
    ax2.scatter(data[:, 0], data[:, 1])
    ax2.plot([5.0, 22.5], [5.0 * w[0] + w[1], 22.5 * w[0] + w[1]], color="#FF0000")
    ax2.set_title("linear regression")
    plt.show()


if __name__ == "__main__":
    maxIter = 10000
    data = np.loadtxt("ex1data1.txt", delimiter=",", dtype=np.float64)

    y = data[:, 1:2]
    X = np.insert(data[:, 0:1], obj=1, values=1, axis=1)
    w = np.zeros(X.shape[1]).reshape(2, 1)
    w, loss = grad_descent(w, X, y, alpha=0.002, maxIter=maxIter)
    print(w)
    paint(loss, data)

最后拟合得到的权重矩阵为$w=\left[\begin{array}{c}1.18218457\\ -3.78778781\end{array}\right]$，损失函数收敛图、线性回归拟合结果如图所示：

多变量线性回归

吴恩达机器学习线性回归例题ex1data2.txt数据集，该数据集共47行×3列，第一列为房子的大小，第二列为卧室数量，第三列为房子的价格，该问题有2个变量（房子的大小，卧室的数量）。请利用线性回归，根据房子的大小和卧室的数量数据预测房子的价格。

2104,3,399900
1600,3,329900
2400,3,369000
1416,2,232000
3000,4,539900
1985,4,299900
1534,3,314900
1427,3,198999
1380,3,212000
1494,3,242500
1940,4,239999
2000,3,347000
1890,3,329999
4478,5,699900
1268,3,259900
2300,4,449900
1320,2,299900
1236,3,199900
2609,4,499998
3031,4,599000
1767,3,252900
1888,2,255000
1604,3,242900
1962,4,259900
3890,3,573900
1100,3,249900
1458,3,464500
2526,3,469000
2200,3,475000
2637,3,299900
1839,2,349900
1000,1,169900
2040,4,314900
3137,3,579900
1811,4,285900
1437,3,249900
1239,3,229900
2132,4,345000
4215,4,549000
2162,4,287000
1664,2,368500
2238,3,329900
2567,4,314000
1200,3,299000
852,2,179900
1852,4,299900
1203,3,239500

数据级 D = { ( x 1 ( i ) , x 2 ( i ) , y ( i ) ) } i = 1 N = 47 \mathcal{D}=\{{(x_1^{(i)},x_2^{(i)}, y^{(i)})} \}_{i=1}^{N=47} D={(x1(i)​,x2(i)​,y(i))}i=1N=47​，$x^{(i)}$为第$i$个样本，即第$i$行的人口数，$x_1$为房子大小，$x_2$为卧室的数量，$N$为样本数量。

加载数据集

导入库：

import numpy as np  			# 科学计算库，处理多维数组，进行数据分析
import matplotlib.pyplot as plt # 提供一个类似 Matlab 的绘图框架

导入数据集：

data = np.loadtxt("./resources/ex1data1.txt", delimiter=',', dtype=np.float64)

特征缩放

在多变量线性回归中，由于变量的范围有所不同，如x1是0到5000的数，而x2则是0到5的数，这会导致在X与Y坐标轴比例尺相同的情况下，整个图像显得极为修长，且此时梯度下降会需要很多次迭代才能收敛。因此需要进行特征缩放，使x1和x2在保持原有对应关系下，能稳定在-1到1之间。
常用特征缩放方法：

1. 线性归一化（Min-Max Normalization）：$X_{norm}=\frac{X-X_{\min }}{X_{\max }-X_{\min }}$把输入数据都转换到[0，1]的范围。
2. 0均值标准化：

假设函数

假设函数为$f(x)=w_1{x}_1+w_2{x}_2+b$，进一步向量乘法的形式，表示起来更简单，易于编程计算，这里构建增广权重向量$w=\left[\begin{array}{c}{w}_1\\ w_2\\ b\end{array}\right]$，增广特征向量$x=\left[\begin{array}{c}{x}_1\\ x_2\\ 1\end{array}\right]$：
KaTeX parse error: {align} can be used only in display mode.

损失函数

令$X=\left[\begin{array}{cc}{x}^{(1)}& 1\\ x^{(2)}& 1\\ ⋮& ⋮\\ x^{(N)}& 1\end{array}\right]$，$y=\left[\begin{array}{c}{y}^{(1)}\\ y^{(2)}\\ ⋮\\ y^{(N)}\end{array}\right]$，则损失函数$\mathcal{L}(w)$可以表示为：

KaTeX parse error: {align} can be used only in display mode.

def loss_func(w, X, y):
    return np.sum(np.square(X.dot(w) - y)) / 2 * X.shape[0]

梯度下降函数

首先对$\mathcal{L}(w)$求导，KaTeX parse error: {align} can be used only in display mode.，这里用到了向量函数及其导数公式：$\frac{\partial Ax}{\partial x}=A^{\mathrm{T}}$。故梯度下降更新公式：
KaTeX parse error: {align} can be used only in display mode.

def grad(w, X, y):
    return (X.T @ (X.dot(w) - y)) / X.shape[0]

def grad_descent(w, X, y, alpha, maxIter):
    loss = []
    for i in range(maxIter):
        w = w - alpha * grad(w, X, y)
        cost = loss_func(w, X, y)
        loss.append(cost)
    return w, loss

算法步骤

1. 初始化最大迭代次数maxIter，学习率alpha；
2. 导入数据：np.loadtxt()，构建假设函数（增广权重矩阵、增广特征向量）$w=\left[\begin{array}{c}0\\ 0\\ 0\end{array}\right]$；
3. 由式（7）计算梯度；
4. 由式（9）更新梯度；
5. 由式（6）计算损失函数；
6. 判断maxIter>当前迭代次数?，是则转step3，否则算法结束，输出增广权重矩阵$w$。

完整代码

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

np.set_printoptions(suppress=True)


def loss_func(w, X, y):
    return np.sum(np.square(X.dot(w) - y)) / 2 * X.shape[0]


def grad(w, X, y):
    return (X.T @ (X.dot(w) - y)) / X.shape[0]


def grad_descent(w, X, y, alpha, maxIter):
    loss = []
    for i in range(maxIter):
        w = w - alpha * grad(w, X, y)
        cost = loss_func(w, X, y)
        loss.append(cost)
    return w, loss


if __name__ == "__main__":
    maxIter = 10000
    # 读取数据
    dataset = np.loadtxt("ex1data2.txt", delimiter=",", dtype=np.float64)
    # 特征缩放
    data = (dataset - np.mean(dataset, axis=0)) / np.std(dataset, axis=0)

    # 增广矩阵
    X = np.insert(data[:, 0:2], obj=2, values=1, axis=1)
    y = data[:, 2:3]

    w = np.zeros(shape=X.shape[1]).reshape(3, 1)
    w, loss = grad_descent(w, X, y, alpha=0.002, maxIter=maxIter)
    print(w)

    fig = plt.figure()
    ax = fig.add_axes(Axes3D(fig))  # 创建三维坐标
    ax.scatter(data[:, 0:1], data[:, 1:2], data[:, -1])  # 散点图
    ax.set(xlabel='size', ylabel='bedrooms', zlabel='price')  # 坐标轴

    # 绘制拟合平面
    x1, x2 = np.meshgrid(
        np.linspace(np.min(data[:, 0:1]), np.max(data[:, 0:1]), 47),
        np.linspace(np.min(data[:, 1:2]), np.max(data[:, 1:2]), 47),
    )  # 生成网格采样点，这一步很重要
    w = w.reshape(1, 3)
    f = w[0][0] * x1 + w[0][1] * x2 + w[0][2]
    ax.plot_surface(x1, x2, f, color='g', alpha=0.4)  # 绘制平面

    plt.show()

运行结果$w=\left[\begin{array}{c}0.88469562\\ -0.05310845\\ -0.\end{array}\right]$，与正规方程w=np.linalg.inv(X.T@X)@X.T@y求解结果一致。