吴恩达机器学习课后作业1.1——多变量线性回归（Linear regression with multiple variable）

1. 问题和数据

假设你要卖掉你的房子，你想知道一个好的市场价格是多少。其中一种方法是，首先收集最近出售的房屋的信息。在本部分的练习中，你将使用多元线性回归来预测房屋价格。
数据ex1data2.txt内容为一个47行3列(47,3)的数据；其中第一列表示房子的面积(单位:平方英尺)，第二列表示卧室的数量，第三列是房子的价格。具体数据如下：

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

2. 数据导入与初步分析

导入包，numpy和pandas是做运算的库，matplotlib是画图的库。

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

导入数据集

path = 'ex1data1.txt'   # 数据路径，要和.py文件放在同一个文件下
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])  #names列名，第一列为Population， 第二列为Profit
print(data.head())  # 预览数据前五行

运行输出结果如下：

      0  1       2
0  2104  3  399900
1  1600  3  329900
2  2400  3  369000
3  1416  2  232000
4  3000  4  539900

此时由于特征值大小不同需要进行均值归一化。
如果这个房子价格不归一化，它的数量级和你输入值规一化数量级差别太大，几十万的数量级和个位小数做回归，就不能保证收敛了，预测的y和实际上y误差几十万差的太多了。

data2 = (data2 - data2.mean()) / data2.std()  # 这是标准通用的均值归一化公式
print(data2.head())

均值归一化之后数据大小就很均匀，均值归一化结果：

          0         1         2
0  0.130010 -0.223675  0.475747
1 -0.504190 -0.223675 -0.084074
2  0.502476 -0.223675  0.228626
3 -0.735723 -1.537767 -0.867025
4  1.257476  1.090417  1.595389

现在使用梯度下降来实现线性回归，以最小化成本函数。
首先，我们将创建一个以参数θ (theta)为特征函数的代价函数。（注：代价函数，成本函数等叫法其实都是一个意思，是$J(\theta )$这个Cost Function。

 

因为该问题是多变量，有两个x，所以本题目中：

 

综合上述两个公式，所以我们建模的公式如下(不会就去看我的“吴恩达机器学习打卡day1——P8”)：

计算代价函数J(θ)

def computeCost(X2, y2, theta):
    inner = np.power(((X2 * theta.T) - y2), 2)  # 即成本函数公式求和符号右边；.T是转置矩阵的意思
    return np.sum(inner) / (2 * len(X2))  # np.sum()可以将矩阵里的元素全部加起来，即成本函数公式求和符号与（1/2m）
# return将它后面的值/结果返回给定义的函数，不然打印函数名出来是空值；有了return之后，所以下面打印函数名才能出现值/结果

让我们在训练集中添加一列，以便我们可以使用向量化的解决方案来计算代价和梯度，
即为了使$h_{\theta }(x)={\theta }_0+{\theta }_{1}X$中${\theta }_0的系数为1$。

#在数据集的最左侧插入一列全为1的列，以便计算
data2.insert(0, 'ones', 1)

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

# .shape会得到该矩阵共有（几行，几列）.shape[0]为第一维的长度，shape[1]为第二维的长度，即列。
# pandas中利用.iloc选取数据loc; ','前的部分标明选取的行，','后的部分标明选取的列， 此时三列了
# set X (training data) and y (target variable)
cols = data2.shape[1]  # 表示得到第二维（行）一共有3列，即cols=3
X2 = data2.iloc[ : , 0:cols-1]  # [X是所有行， y是从0到cols-1列（即第一、二列）]
y2 = data2.iloc[ : , cols-1: cols]  # [X是所有行， y是从cols-1到cols列（即第三列）]

打印X表头，观察下X（训练集）是否获取正确。

print(X2.head())

输出X表头：

   ones         0         1
0     1  0.130010 -0.223675
1     1 -0.504190 -0.223675
2     1  0.502476 -0.223675
3     1 -0.735723 -1.537767
4     1  1.257476  1.090417

打印y表头，观察下y（训练集）是否获取正确。

print(y2.head())
```python

输出y表头：

          2
0  0.475747
1 -0.084074
2  0.228626
3 -0.867025
4  1.595389

代价函数是应该是numpy矩阵，所以我们需要转换x和y成矩阵，然后才能使用它们。 我们还需要初始化theta，即把theta所有元素都设置为0。

# pandas读取的数据是DataFrame形式，优势是能对数据进行很多操作,但要想进行矩阵运算，要将DataFrame形式转换为矩阵，例如x=np.matrix(X.values)
X2 = np.matrix(X2.values)  # 将X转换为矩阵
y2 = np.matrix(y2.values)  # 将y转换为矩阵
# theta = np.matrix(np.array([0, 0, 0]))  # 因为X变成了三列，所以theta自然也要有三列
theta = np.full((1, X2.shape[1]), 0)   # np.full()是填充0的函数，与np.zeros()差不多; 79行的.ravel去掉了我觉得用它没啥用
# np.full((行数，列数)， 想要填充的数字)

#打印theta检查，应该是个一行两列的零矩阵

print(theta, 'theta')

输出theta：

theta: [[0 0 0]]

看下维度，维度对了矩阵的乘法计算才能正确进行，这一作业的代码错误经常都是矩阵维度不对造成的。

print('X.shape:', X.shape)
print('theta.shape', theta.shape)
print('y.shape', y.shape)

输出X2, theta, y2的维度：

X2.shape: (47, 3)
theta.shape (1, 3)
y2.shape (47, 1)

计算初始的代价函数（theta初始值为0）

computeCost(X, y, theta)
print('Cost_init:', computeCost(X, y, theta))  #得到初始的成本代价（还未开始迭代）

初始代价函数结果

Cost_init: 32.072733877455676

3. 批量梯度下降（Batch gradient decent）

关键点在于theta0和thata1要同时更新，需要用temp进行临时存储。
temp[0, j]这一行是$\theta$的迭代公式，也是梯度下降函数的核心。也就是：

def gradientDescent(X2, y2, theta, alpha, iters):
    temp = np.matrix(np.zeros(theta.shape))   # 按照theta的shape建立一个（1，3）的零值矩阵temp，为了临时存储theta，方便进行迭代更新。
    parameters = int(theta.shape[1])  # theta.shape[1]就是3; ravel计算需要求解的参数个数，功能将多维数组降至一维
    print(parameters, 'parameters')
    cost = np.zeros(iters) # 构建iters个0的数组，用来存放cost

    for i in range(iters):
        error = (X2 * theta.T) - y2  # 计算出每一组数据的误差，每行都是样本的误差

        for j in range(parameters):
            # 线性回归方程对theta求导后的公式为(theta.T*X-y)*xj的平方的sum
            # 所以这里用X*theta.T，即列向保存，再乘以X中对应的每一列，方便保存求和。
            # multiply为数组对应元素相乘，所以这个term就是还未相加的对theta求导后的数据
            term = np.multiply(error, X2[ : , j])
            temp[0, j] = theta[0, j] - ((alpha / len(X2)) * np.sum(term))  # 前面转不转置和上课学的公式不太一样，这没关系，因为书
            # 上的转置是考虑了求和，我们这里的转置是为了变成可以相乘的方块，至于求和最后用了sum直接求

        theta = temp
        cost[i] = computeCost(X2, y2, theta)  # 每迭代一次就调用一次成本函数，计算目标方程的值；并存为cost数组里面的第i个值

    return theta, cost

初始化学习速率和迭代次数

alpha = 0.01
iters = 1000

#调用梯度下降函数，计算迭代后的最优值

g2, cost2 = gradientDescent(X2, y2, theta, alpha, iters)
# g为迭代完之后的theta;（ 按return的theta， cost排序）
print(g2, "g2")  # 此时的g为满足使成本函数最小的最优值theta
print(cost2, "cost2")

输出最优g2为：

[[-1.10910099e-16  8.78503652e-01 -4.69166570e-02]] g2

代入最优值theta，计算最小成本

minCost = computeCost(X2, y2, g2)  # 代入最优值theta，计算最小成本函数
print(minCost, "minCost")

输出最小成本为：

0.1307033696077189 minCost

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

# fig代表绘图窗口(Figure)；ax代表这个绘图窗口上的坐标系(axis)，一般会继续对ax进行操作
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(np.arange(iters), cost2, 'r')  # np.arange()自动返回等差数组
ax.set_xlabel('Iterations')  # 设定x轴名称
ax.set_ylabel('Cost')  # 设定y轴名称
ax.set_title('Error vs. Training Epoch')  # 设定整个表的名称
plt.show()

输出绘制的迭代曲线：

参考文献：
[1] https://blog.csdn.net/weixin_43455338/article/details/104794760
[2] https://github.com/PlayPurEo/ML-and-DL/tree/master/basic-model
[3] https://www.bilibili.com/video/BV1cX4y1G7h2?spm_id_from=333.999.0.0
[4] https://www.bilibili.com/video/BV1Xt411s7KY?spm_id_from=333.999.0.0