1. 问题和数据
假设你是一家连锁餐车店的老板，但是你又和别的土老板不一样，你又刚好是个懂线性回归，还懂编程的老板，正在考虑在不同的城市开一家新店。该连锁店已经在各个城市开设了餐车，你可以获得这些城市的利润和人口数据。在本部分的练习中，您将透过这些数据，通过单变量实现线性回归，以预测餐车的利润。
数据ex1data1.txt内容为一个97行两列(97,2)的数据；其中第一列表示人口数据，第二列表示利润。具体数据如下：
6.1101,17.5925.5277,9.13028.5186,13.6627.0032,11.8545.8598,6.82338.3829,11.8867.4764,4.34838.5781,126.4862,6.59875.0546,3.81665.7107,3.252214.164,15.5055.734,3.15518.4084,7.22585.6407,0.716185.3794,3.51296.3654,5.30485.1301,0.560776.4296,3.65187.0708,5.38936.1891,3.138620.27,21.7675.4901,4.2636.3261,5.18755.5649,3.082518.945,22.63812.828,13.50110.957,7.046713.176,14.69222.203,24.1475.2524,-1.226.5894,5.99669.2482,12.1345.8918,1.84958.2111,6.54267.9334,4.56238.0959,4.11645.6063,3.392812.836,10.1176.3534,5.49745.4069,0.556576.8825,3.911511.708,5.38545.7737,2.44067.8247,6.73187.0931,1.04635.0702,5.13375.8014,1.84411.7,8.00435.5416,1.01797.5402,6.75045.3077,1.83967.4239,4.28857.6031,4.99816.3328,1.42336.3589,-1.42116.2742,2.47565.6397,4.60429.3102,3.96249.4536,5.41418.8254,5.16945.1793,-0.7427921.279,17.92914.908,12.05418.959,17.0547.2182,4.88528.2951,5.744210.236,7.77545.4994,1.017320.341,20.99210.136,6.67997.3345,4.02596.0062,1.27847.2259,3.34115.0269,-2.68076.5479,0.296787.5386,3.88455.0365,5.701410.274,6.75265.1077,2.05765.7292,0.479535.1884,0.204216.3557,0.678619.7687,7.54356.5159,5.34368.5172,4.24159.1802,6.79816.002,0.926955.5204,0.1525.0594,2.82145.7077,1.84517.6366,4.29595.8707,7.20295.3054,1.98698.2934,0.1445413.394,9.05515.4369,0.61705
2. 数据导入与初步分析
导入包，numpy和pandas是做运算的库，matplotlib是画图的库。
import numpy as npimport pandas as pdimport matplotlib.pyplot as plt
导入数据集
path = 'ex1data1.txt' # 数据路径，要和.py文件放在同一个文件下data = pd.read_csv(path, header=None, names=['Population', 'Profit']) #names列名，第一列为Population， 第二列为Profitprint(data.head()) # 预览数据前五行
运行输出结果如下：
Population Profit06.1101 17.59.5277 9.130228.5186 13.662037.0032 11.854045.8598 6.8233
可以描述一下数据值。
对于数值数据，描述的结果将包括计数，平均值，标准差，最小值，最大值以及较低的百分位数和50。默认情况下，较低的百分位数为25，较高的百分位数为75，百分位数与中位数相同。
print(data.describe())
描述结果：
PopulationProfitcount 97.000000 97.000000mean8.159800 5.839135std3.869884 5.510262min5.026900 -2.68070025%5.707700 1.98690050%6.589400 4.56230075%8.578100 7.046700max22.203000 24.147000
在开始任何任务之前，通过可视化来理解数据通常是有用的。
对于这个数据集，您可以使用散点图来可视化数据，因为它只有两个属性(利润和人口)。
# 数据可视化，kind选择绘制散点图，取值为line或者scatter；依次选择x列的名字，y列的名字；figsize默认值图像大小data.plot(kind='scatter', x='Population', y='Profit', figsize=(12, 8))plt.show()
输出的散点图如下，可以直接看出x和y还是有着较强的线性关系的，这就是画出散点图的好处。
现在使用梯度下降来实现线性回归，以最小化成本函数。
首先，我们将创建一个以参数θ (theta)为特征函数的代价函数。（注：代价函数，成本函数等叫法其实都是一个意思，是J(θ)J(\theta)J(θ)这个Cost Function。
因为该问题是单变量，只有一个x，所以本题目中：
综合上述两个公式，所以我们建模的公式如下(不会就去看我的“吴恩达机器学习打卡day1——P8”)：
计算代价函数J(θ)
def computeCost(X, y, theta):inner = np.power(((X * theta.T) - y), 2) # 即代价函数公式求和符号右边；.T是转置矩阵的意思return np.sum(inner) / (2 * len(X)) # np.sum()可以将矩阵里的元素全部加起来，(2 * len(X))是成本函数公式求和符号与（1/2m）# return将它后面的值/结果返回给定义的函数，不然打印函数名出来是空值；有了return之后，所以下面打印函数名才能出现值/结果
让我们在训练集中添加一列，以便我们可以使用向量化的解决方案来计算代价和梯度，
即为了使hθ(x)=θ0+θ1Xh_\theta(x) = \theta_0 +\theta_1 Xhθ(x)=θ0+θ1X中θ0的系数为1\theta_0的系数为1θ0的系数为1。
#在数据集的最左侧插入一列全为1的列，以便计算data.insert(0, 'ones', 1) # 即x0=1 索引index为O, name为ones,value为1.
现在我们来做一些变量初始化。
# .shape会得到该矩阵共有（几行，几列）.shape[0]为第一维的长度，shape[1]为第二维的长度，即列。# pandas中利用.iloc选取数据loc; ','前的部分标明选取的行，','后的部分标明选取的列， 此时三列了# set X (training data) and y (target variable)cols = data.shape[1] # 表示得到第二维（行）一共有3列，即cols=3X = data.iloc[ : , 0:cols-1] # [X是所有行， y是从0到cols-1列（即第一、二列）]y = data.iloc[ : , cols-1: cols] # [X是所有行， y是从cols-1到cols列（即第三列）]
打印X表头，观察下X（训练集）是否获取正确。
print(X2.head())
输出X表头：
ones Population016.1101115.5277218.5186317.0032415.8598
打印y表头，观察下y（训练集）是否获取正确。
print(y2.head())```python
输出y表头：
Profit0 17.59201 9.13022 13.66203 11.85404 6.8233
代价函数是应该是numpy矩阵，所以我们需要转换x和y成矩阵，然后才能使用它们。 我们还需要初始化theta，即把theta所有元素都设置为0。
# pandas读取的数据是DataFrame形式，优势是能对数据进行很多操作,但要想进行矩阵运算，要将DataFrame形式转换为矩阵，例如x=np.matrix(X.values)X = np.matrix(X.values) # 将X转换为矩阵y = np.matrix(y.values) # 将y转换为矩阵theta = np.matrix(np.array([0, 0])) # 因为X是两列，所以theta自然也要有两列
#打印theta检查，应该是个一行两列的零矩阵
print(theta，'theta')
输出theta：
theta: [[0 0]]
看下维度，维度对了矩阵的乘法计算才能正确进行，这一作业的代码错误经常都是矩阵维度不对造成的。
print('X.shape:', X.shape)print('theta.shape', theta.shape)print('y.shape', y.shape)
计算初始的代价函数（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(X, y, theta, alpha, iters):temp = np.matrix(np.zeros(theta.shape)) # 按照theta的shape建立一个（1，2）的零值矩阵temp，为了临时存储theta，方便进行迭代更新。parameters = int(theta.shape[1]) # theta.shape[1]就是2，parameter就是$\theta_j$的j的意思;cost = np.zeros(iters) # 构建iters个0的数组，用来存放costfor i in range(iters):error = (X * theta.T) - y # 计算出每一组数据的误差，每行都是样本的误差for j in range(parameters):# 线性回归方程对theta求导后的公式为(theta.T*X-y)*xj的平方的sum# 所以这里用X*theta.T，即列向保存，再乘以X中对应的每一列，方便保存求和。# multiply为数组对应元素相乘，所以这个term就是还未相加的对theta求导后的数据term = np.multiply(error, X[ : , j]) # X[ : , j]的意思是X[所有行，第j列]temp[0, j] = theta[0, j] - ((alpha / len(X)) * np.sum(term)) # 前面转不转置和上课学的公式不太一样，这没关系，因为书# 上的转置是考虑了求和，我们这里的转置是为了变成可以相乘的方块，至于求和最后用了sum直接求theta = tempcost[i] = computeCost(X, y, theta) # 每迭代一次就调用一次成本函数，计算目标方程的值；并存为cost数组里面的第i个值return theta, cost
初始化学习速率和迭代次数
alpha = 0.01 # 可以试试0.005，0.001越小迭代得越快iters = 1000 # 10000时，迭代也曲线挺好的
#调用梯度下降函数，计算迭代后的最优值
g, cost = gradientDescent(X, y, theta, alpha, iters)# g为迭代完之后的theta;（ 按return的theta， cost排序）print('g:', g) # 此时的g为满足使成本函数最小的最优值theta
输出最优g为：
g： [[-3.24140214 1.1272942 ]]
代入最优值theta，计算最小成本
minCost = computeCost(X, y, g) # 代入最优值theta，计算最小成本函数print('minCost:', minCost)
输出最小成本为：
minCost: 4.515955503078914
现在我们来绘制纯属模型以及数据，以便直观地看出它的拟合
x = np.linspace(data.Population.min(), data.Population.max(), 100) # 从最小数据到最大数据（即所有数据）中随机抽100个数据点，# 用np.linspace可以将这些点用等差数列排列出来，即得到一个个间隔均匀的x点# 线性回归的最终方程f = g[0, 0] + (g[0, 1] * x) # 将一个个间隔均匀的x点带入h(\theta)函数，得到一个个f值,也就是我们的预测值fig, ax = plt.subplots(figsize=(12, 8)) # fig代表绘图窗口(Figure)；ax代表这个绘图窗口上的坐标系(axis)，一般会继续对ax进行操作。figsize用来设置图形的大小ax.plot(x, f, 'r', label='Prediction') # 设置横纵坐标的函数，并设置颜色为红色，在图标上标上'Prediction'标签ax.scatter(data.Population, data.Profit, label='Training Data') # 画散点图，设置横纵坐标的函数，并设置颜色为红色，在图上标上'Training Data'标签ax.legend(loc=4) # 给标签选个位置，1表示在第一象限（右上角），2表示在左上角，3表示在左下角，4表示在第四象限ax.set_xlabel('Population') # 设定x轴名称ax.set_ylabel('Profit') # 设定y轴名称ax.set_title('Predicted Profit vs. Population Size') # 设定整个表的名称plt.show()
输出绘制的拟合曲线为：.
由于梯度方程函数在每个训练迭代中输出一个代价的向量，所以我们也可以绘制代价函数的迭代曲线。请注意，代价总是在降低的，这是凸优化问题的一个特点。
代价函数的曲线
#fig代表绘图窗口(Figure)；ax代表这个绘图窗口上的坐标系(axis)，一般会继续对ax进行操作fig, ax = plt.subplots(figsize=(12, 8))ax.plot(np.arange(iters), cost, 'r') # np.arange()自动返回等差数组ax.set_xlabel('Iterations') # 设定x轴名称ax.set_ylabel('Cost') # 设定y轴名称ax.set_title('Error vs. Training Epoch') # 设定整个表的名称plt.show()
输出绘制的迭代曲线：
将迭代次数iters从1000换成10000，迭代更快，迭代曲线更优一点：
参考文献：
[1] /weixin_43455338/article/details/104794760
[2] /PlayPurEo/ML-and-DL/tree/master/basic-model
[3] /video/BV1cX4y1G7h2?spm_id_from=333.999.0.0
[4] /video/BV1Xt411s7KY?spm_id_from=333.999.0.0
