练习一 - 线性回归(Linear regression)

一、引言

在本次练习中，我们将会学习到线性回归和及其工作模式。当然，本次练习是基于课程内容的，强烈建议看完课程后再完成本章内容。

另外，本次练习需要调用到数据分析的一些包：Numpy，Matplotlib， Pandas等；需要将这些数据分析包的基本操作了解一下(菜鸟教程)。同时我们需要将所用数据(ex1data1.txt,ex1data2.txt)导入进来，可以去GitHub上获取。

二、单变量线性回归

1. 问题背景

本部分中，我们将会实现一个可以用于预测副食卡车利润的单变量线性回归模型。假如你是一个餐店零售商的CEO，正在考虑在不同的城市中增设一些新档口；并且你已拥有目前各个城市中不同连锁店的利润和人流数据。

文件ex1data1.txt包含了解决本次线性回归问题的数据集。第一列是某个城市连锁店的人流量，第二列是该城市的连锁店的利润(负值则为亏损)。

2. 导入数据

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd

# 使用pandas导入数据，将数据集使用dataFrame结构导入
# 数据集路径
path =  'ex1data1.txt'
# header: 指定第几行作为列名(忽略注解行)，如果没有指定列名，默认header=0; 如果指定了列名header=None
# names 指定列名，如果文件中不包含header的行，应该显性表示header=None
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
data

# pandas中的DataFrame相当于二维数组
type(data)

pandas.core.frame.DataFrame

3. 绘制数据

data.describe()

• DataFrame.describe(percentiles=None,include=None,exclude=None)

1. 作用：用于生成描述性统计数据，统计数据集的集中趋势，分散和行列的分布情况，不包括 NaN值。
2. 参数：

• percentiles：赋值类似列表形式，可选表示百分位数，介于0和1之间。默认值为 [.25,.5,.75]，分别返回第25，第50和第75百分位数。 可自定义其它值，用法为df.describe(percentiles=[.xx])。

• include：'all'，类似于dtypes列表或None（默认值），可选要包含在结果中的数据类型的白名单。对于Series不可用。以下是选项：'all'：输入的所有列都将包含在输出中。类似于dtypes的列表：将结果限制为提供的数据类型。将结果限制为数字类型用法：numpy.number。要将其限制为对象列用法：numpy.object。字符串也可以以select_dtypes（例如df.describe(include=['O'])）的方式使用。要选择分类类型，请使用'category'无（默认）：结果将包括所有数字列。

• exclude：类似于dtypes列表或None（默认值），可选，要从结果中除去的黑名单数据类型列表。Series不可用。以下是选项：类似于dtypes的列表：从结果中排除提供的数据类型。排除数值类型用法：numpy.number。要排除对象列，使用numpy.object。字符串也可以以select_dtypes（例如df.describe(include=['O'])）的方式使用。要排除分类类型，请使用'category'无（默认）：结果将不包含任何内容。

3. 注意：percentiles只会在数据有序排列时才会有分析作用

# 使用matplotlib + pandas绘制数据
# 参数含义分别为：散点图、x轴、y轴，画布大小
data.plot(kind = 'scatter', x = 'Population', y = 'Profit',marker = 'x' ,color = 'r',figsize = (12., 8))
plt.show()

三、批量梯度下降

现在让我们使用梯度下降来实现线性回归，以最小化成本函数。 以下代码示例中实现的方程在“练习”文件夹中的“ex1.pdf”中有详细说明。

首先，我们将创建一个以参数θ为特征函数的代价函数

$J(\theta) = \frac{1}{2m}\sum^m_{ i =1}(h_\theta(x^{(i)}) - y^{(i)})^2 \\ h_\theta(x) = \theta^T X = \theta_0 x_0 + \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n$

此处的m毋庸置疑就是data的行数 —— 97,n就是特征向量的长度也就是x的列数 -1 ,因为最后一列是标签y

# my method
def computeCost(X, y, theta):
    res = 0
    # 此处的 i 用于遍历每一行，上面的X就是第(i)行的特征向量
    for i in range(len(X)):
        res += np.power(((X[i] * theta.T) - y[i]), 2)
    return res.sum() / (2 * len(X))
# 简洁表示
# def computeCost(X, y, theta):
#     inner = np.power(((X * theta.T) - y), 2)
#     return np.sum(inner) / (2 * len(X))

回想一下，您的模型的参数是 $\theta_j$值。 这些$\theta$是我们需要优化来使代价函数$J(\theta)$值最小化的值。一种解决办法就是使用批量梯度下降算法。 在批量梯度下降中，每个迭代执行如下更新

$\theta_{j}:=\theta_{j}-\alpha \frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}$

我们是同时更新所有$\theta_j$值的，也就是更新参数向量$\theta$

• 注意：我们将样本数据集逐个存储在X阵列中的每一行了，但是为了便于实现(主要是向量乘法)我们将参数$\theta_0$列为考虑，所以我们在X中添加了额外一列$(x_0^{(i)})$并将其置为全一，那么我们在实现这个模型时，就可将$\theta_0^{(i)}$视为另一个"特征"

1. 数据初始化

之前我们已经使用pandas + matplotlib 导入数据并绘制，接下来我们需要在data中加入一列数据并置为全一,以便我们可以使用向量化的解决方案来计算代价和梯度。

# 插入列，表头，值
data.insert(0, 'Ones', 1)

# head()是观察前5行
data.head()

现在我们来做一些变量初始化。

# 获取列数
cols = data.shape[1]
# X取自data中所有行，去掉最后一列 
X = data.iloc[:,0:cols-1]
# y取自data中所有行，最后一列
y = data.iloc[:,cols-1:cols]

X.head()

y.head()

代价函数是应该是numpy矩阵，所以我们需要转换X和Y(这是pandas的DataFrame数据)，然后才能使用它们。 我们还需要初始化theta。

type(X)

pandas.core.frame.DataFrame

X = np.matrix(X.values)
y = np.matrix(y.values)
theta = np.matrix(np.zeros(X.shape[1]))

theta 是一个(1,2)矩阵

theta

matrix([[0., 0.]])

观察一下各参数的维度

X.shape, theta.shape, y.shape

((97, 2), (1, 2), (97, 1))

首次计算代价函数 (theta初始值为0)。那么代价函数计算的就是：$J(\theta) = \frac{1}{2m} \sum^m_{i = 1}( 0 - y^{(i)} )^2$

# 注意数据转换，否则res为向量而非标量
computeCost(X, y, theta)

32.072733877455654

# 可以验证
res = np.power(y,2)
res = np.sum(res)/ (2 * len(y))
res

32.072733877455676

2. 批量梯度下降

之前我们分析过，我们的最终目的就是得到是使代价函数$J(\theta)$最小的向量$\theta$，由于$J(\theta)$是一个凸函数，那么我们在每次梯度下降时，我们都会进行如下更新：

$\theta_{j}:=\theta_{j}-\alpha \frac{1}{m} \sum_{i=1}^{m}\left(h_{\theta}\left(x^{(i)}\right)-y^{(i)}\right) x_{j}^{(i)}$

定义梯度下降函数

# 参数分别为：训练集及其标签，参数(权重)，学习率(步长)，【超参数】迭代次数
def gradientDescent(X, y, theta, alpha, iters):
    # 定义工具θ矩阵，用于存储中间结果
    temp = np.matrix(np.zeros(theta.shape))
    # 获取参数向量 Θ 的大小
    parameters = int(theta.ravel().shape[1])
    # 用于存储每次迭代的代价函数结果的一维阵列
    cost = np.zeros(iters)
    
    for i in range(iters):
        # 计算误差阵列，error的维度和y保持一致
        error = (X * theta.T) - y
        # 同时更新所有θ值，也就是整个参数矩阵，j的取值范围为 [0, 1],两列
        for j in range(parameters):
            term = np.multiply(error, X[:,j])
            temp[0,j] = theta[0,j] - ((alpha / len(X)) * np.sum(term))
        # 迭代参数矩阵 θ并记录此次迭代的代价函数值    
        theta = temp
        cost[i] = computeCost(X, y, theta)
        
    return theta, cost

初始化一些附加变量 - 学习速率α和要执行的迭代次数。

alpha = 0.01
iters = 1000

现在让我们运行梯度下降算法来将我们的参数θ适合于训练集，并使代价函数最小化。

g, cost = gradientDescent(X, y, theta, alpha, iters)
# 查看最终参数矩阵
g

matrix([[-3.24140214,  1.1272942 ]])

最后，我们可以使用我们拟合的参数计算训练模型的代价函数（误差）。

computeCost(X, y, g)

4.515955503078912

可以看出，我们使用梯度下降得到的θ确实使代价函数值减小不少(与之前摆烂的 0 值相比)

现在我们来绘制线性模型以及数据，直观地看出它的拟合。

3. 绘制模型

# 正规化x轴
x = np.linspace(data.Population.min(), data.Population.max(), 100)
# 得到假设函数直线
f = g[0, 0] + (g[0, 1] * x)

# 预备组合图
fig, ax = plt.subplots(figsize=(12,8))
# 绘制最终假设函数——模型直线
ax.plot(x, f, 'r', label='Prediction')
# 绘制样本数据集
ax.scatter(data.Population, data.Profit, label='Traning Data')

ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()

由于梯度方程式函数也在每个训练迭代中输出一个代价的向量，所以我们也可以绘制。 请注意，代价一定是逐步降低的 - 这是凸优化问题的一个例子。

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

四、多变量线性回归

练习1还包括一个房屋价格数据集(ex1data2.txt)，其中有2个变量（房子的大小，卧室的数量）和目标（房子的价格）。 我们使用我们已经应用的技术来分析数据集。

1. 导入数据

path =  'ex1data2.txt'
data2 = pd.read_csv(path, header=None, names=['Size', 'Bedrooms', 'Price'])
data2.head()

对于此任务，我们添加了另一个预处理步骤 —— 特征归一化。 这个对于pandas来说很简单

data2 = (data2 - data2.mean()) / data2.std()
data2.head()

2. 数据预处理并训练

现在我们重复第1部分的预处理步骤，并对新数据集运行线性回归程序。

# 添加额外一列并置一
data2.insert(0, 'Ones', 1)

# 分离训练样本中的特征和标签
cols = data2.shape[1]
X2 = data2.iloc[:,0:cols-1]
y2 = data2.iloc[:,cols-1:cols]

# 将pandas数据转换为numpy矩阵
X2 = np.matrix(X2.values)
y2 = np.matrix(y2.values)
# 由于X的特征向量(未填充)的长度为 2，所以参数向量的长度为 3
theta2 = np.matrix(np.array([0,0,0]))

# 基于训练集训练线性回归
g2, cost2 = gradientDescent(X2, y2, theta2, alpha, iters)

# 计算最终模型的代价函数
computeCost(X2, y2, g2)

0.13070336960771892

我们也可以快速查看这一个的训练进程。

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost2, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

我们也可以使用scikit-learn的线性回归函数，而不是从头开始实现这些算法。 我们将scikit-learn的线性回归算法应用于第1部分的数据，并看看它的表现。

五、scikit-learn训练线性回归

from sklearn import linear_model
model = linear_model.LinearRegression()
# 可能会报警告，是由于numpy本本过高引起的，后面报错可以降一下版本
model.fit(X, y)

E:\anaconda\lib\site-packages\sklearn\utils\validation.py:593: FutureWarning: np.matrix usage is deprecated in 1.0 and will raise a TypeError in 1.2. Please convert to a numpy array with np.asarray. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.matrix.html
  warnings.warn(
E:\anaconda\lib\site-packages\sklearn\utils\validation.py:593: FutureWarning: np.matrix usage is deprecated in 1.0 and will raise a TypeError in 1.2. Please convert to a numpy array with np.asarray. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.matrix.html
  warnings.warn(





LinearRegression()

scikit-learn model的预测表现

# 获取所有样本(行)的特征值(人流量——第二列)
x = np.array(X[:, 1].A1)
# 得到假设函数
f = model.predict(X).flatten()

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x, f, 'r', label='Prediction')
ax.scatter(data.Population, data.Profit, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()

E:\anaconda\lib\site-packages\sklearn\utils\validation.py:593: FutureWarning: np.matrix usage is deprecated in 1.0 and will raise a TypeError in 1.2. Please convert to a numpy array with np.asarray. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.matrix.html
  warnings.warn(

六、normal equation（正规方程）

1. 方法

正规方程是通过求解下面的方程来找出使得代价函数最小的参数的： $\frac{\partial}{\partial \theta_j}J(\theta_j) = 0$

假设我们的训练集特征矩阵为 X(填充了$x_0^{(i)} = 1$)之后，并且我们的训练集标签为向量y，则我们利用正规方程解出的向量 $\theta$ 即为所求： $\theta = (X^TX)^{-1}X^Ty$

其中， 上标T代表矩阵转置，上标-1 代表矩阵的逆。

设矩阵$A = X^TX$，则：$(X^TX) ^{-1} = A^{-1}$

2. 梯度下降与正规方程的比较

1. 梯度下降：需要选择学习率α，需要多次迭代，当特征数量n大时也能较好适用，适用于各种类型的模型

2. 正规方程：不需要选择学习率α，一次计算得出，需要计算 $(𝑋^𝑇𝑋)^{−1}$，如果特征数量n较大则运算代价大，因为矩阵逆的计算时间复杂度为 𝑂(𝑛^3) ，通常来说当 𝑛 小于10000 时还是可以接受的，只适用于线性模型，不适合逻辑回归模型等其他模型

3. 运用

# 正规方程
def normalEqn(X, y):
    #X.T@X等价于X.T.dot(X) ，矩阵乘法
    theta = np.linalg.inv(X.T@X)@X.T@y
    # 返回的是参数向量
    return theta

eq_theta2=normalEqn(X, y).reshape(1,2)
eq_theta2

matrix([[-3.89578088,  1.19303364]])

# 梯度下降得到的结果
ga_theta,cost = gradientDescent(X, y, theta, alpha, iters)
ga_theta

matrix([[-3.24140214,  1.1272942 ]])

# 最终代价函数
# 注意转变参数矩阵的尺寸
eq_finalcost = computeCost(X,y,eq_theta2)
ga_finalcost = computeCost(X,y,ga_theta)
print('梯度下降得到的最终代价函数结果：%f \n正规方程得到的最终代价函数结果：%f'%(ga_finalcost, eq_finalcost))

梯度下降得到的最终代价函数结果：4.515956 
正规方程得到的最终代价函数结果：4.476971

这样看来，正规方程在数据规模比较小的时候计算得又快又准！

在练习2中，我们学习如何解决分类问题的逻辑回归。