机器学习代码实现：多元线性回归（梯度下降法）吴恩达课后题目

数据集

数据集如图：（面积，卧室数，价格），来自机器学习吴恩达的课后作业

预测价格

数据集下载：链接: /s/1MzUq1jPVlic5kkTGsXY87Q?pwd=hdkk 提取码: hdkk

–来自百度网盘超级会员v4的分享

代码思路

1、模型思想

当样本xi⃗=(xi1,xi2,..,xid)\vec{x_i} = (x_{i1}, x_{i2},..,x_{id})xi=(xi1,xi2,..,xid)有多个属性描述时，我们采用多元线性回归模型，使得模型预测值f(xi⃗)=w⃗Txi⃗+bf(\vec{x_i}) = \vec{w}^T\vec{x_i} + bf(xi)=wTxi+b，与真实标记yiy_iyi之间的差距尽可能小。

代码中XXX如图：即f(xi⃗)=w⃗^Txi⃗^f(\vec{x_i}) = \hat{\vec{w}}^T \hat{\vec{x_i}}f(xi)=w^Txi^，给每个xi⃗\vec{x_i}xi添加一列1，省略b的计算。

具体公式推导理论这里不多说，详情看西瓜书，之后会写一篇总结。

2、策略：代价函数的选取

依然是差距平方和最小化：

此计算代价的函数不管一元还是多元都是不变的：

# 代价函数，返回该模型w，b参数下的代价def cost_function(w_hat_matrix, X_matrix, y_matrix) -> float:""" 此方法不管几维都固定不变:param w_hat_matrix: 二维matrix，1 * 3:param X_matrix: 二维matrix，47 * 3:param y_matrix: 二维matrix，47 * 1:return: 代价"""m = len(X_matrix) # 样本数return np.sum(np.power(X_matrix * w_hat_matrix.T - y_matrix, 2)) / (2 * m)

3、计算w和b：梯度下降法

本题采用梯度下降法实现，可以得到局部最优解。

梯度下降法可以参考吴恩达讲解的，暂时会个公式就可以：注意对代价函数求偏导

1、吴恩达讲的θ\thetaθ就是西瓜书里的www，h就是fh就是fh就是f

2、每梯度下降一次，会使得代价更小

3、影响每次下降的因素有：学习率α\alphaα和初始的www。

一般初始w⃗^=[[0,0,...应d+1个0]]\hat{\vec{w}} = [[0,0,...应d+1个0]]w^=[[0,0,...应d+1个0]]，即1行d+1列，d为属性个数学习率α\alphaα 一般设为0.01，或者0.03， 0.1， 0.3， 1， 3， 10。

4、迭代指定次数后，返回：

最后一次迭代得到的w⃗^\hat{\vec{w}}w^,将此结果带入代价函数中，计算得到最小的costcost一维array，初始化cost = np.zeros(iters),记录每次迭代的代价，最后画出曲线，查看收敛速度。

代码：梯度下降

# 梯度下降函数def gradientDescent(alpha, iters, w_hat_matrix, X_matrix, y_matrix) -> tuple:""":param alpha: 梯度下降学习率/步长:param iters: 梯度下降次数:param w_hat_matrix: 二维matrix，1 * 3；一般初始设置为[[0, 0, 0]]，初值对梯度下降收敛速度影响大:param X_matrix: 二维matrix，47 * 3。就是公式推导里记的那个超大矩阵X.存放扩展了1的所有样本:param y_matrix: 二维matrix，47 * 1。存放真实标记:return:"""parameters = int(w_hat_matrix.shape[1]) # 3个参数m = len(X_matrix) # 样本数47cur_w_mat_matrix = np.matrix(np.zeros(w_hat_matrix.shape)) # 暂存每次迭代得到的w_hat1*3cost = np.zeros(iters) # 记录每次迭代后的新的代价cost# 迭代iters次for i in range(iters):# 所有样本x_i预测输出和真实标记 y_i的误差。第i行记录了样本x_i的误差error_matrix = X_matrix * w_hat_matrix.T - y_matrix # 97*1# 梯度下降公式：得到新的w_hat：j指向每一列，进行更新（w_1, w_2,..,w_j,...w_d, b）for j in range(parameters):# np.multiply：m*n和m*n的相同下标元素元素相乘，结果还是m*n的矩阵。可以看成多个数乘s = np.sum(np.multiply(error_matrix, X_matrix[:, j])) # 公式：error矩阵和x的第j列相乘cur_w_mat_matrix[0, j] = cur_w_mat_matrix[0, j] - alpha * s / mw_hat_matrix = cur_w_mat_matrix # 每次下降更新w_hatcost[i] = cost_function(w_hat_matrix, X_matrix, y_matrix) # 更新代价，用于观察梯度下降的代价变化曲线return w_hat_matrix, cost

4、特征缩放

4.1 均值标准化

使用均值标准化，因为样本中两个属性值相差太大，代价函数收敛会很慢，因此要先对数据进行缩放。

采用公式：data−平均值标准差data - 平均值 \over 标准差标准差data−平均值

data = pd.read_csv(path, names=["Sizes", "Bedrooms", "Prices"])# 获取原始数据data的一些描述data_origin = data.values # 二维array，每个子array存放样例的数据means = data.mean().values # 一维array，第i个元素：第i列的均值stds = data.std().values # 一维array，第i个元素：第i列的标准差mins = data.min().values # 一维array，第i个元素：第i列的最小值maxs = data.max().values # 一维array，第i个元素：第i列的最大值'''data特征缩放：均值标准化，让不同特征值差异不要太大，否则梯度下降收敛会很慢'''data = (data - data.mean()) / data.std()

4.2 恢复

"""之前获得的res_w_hat，是数据缩放后得到的res_w_hat最终我们要把res_w_hat恢复成与原始数据对应的transform_w_hatdata = (data - data.mean()) / data.std()"""# 咋说呢，就是做了个转置res_w_hat_reshape = np.array(res_w_hat.reshape(-1, 1)) # 把1*3的二维matrix变成 3*1 的二维arraymeans_reshape = means.reshape(-1, 1) # 3 * 1 的二维arraystds_reshape = stds.reshape(-1, 1) # 3*1 的二维arraytransform_w_hat = w_hat_transform(res_w_hat_reshape, means_reshape, stds_reshape)print("transform_w_hat:", transform_w_hat)

# 把数据经过特征缩放（均值标准化）的w_hat变成符合原始数据的w_hat_transformdef w_hat_transform(arr_w_hat_T, means_T, stds_T):""" 同型array相乘相除：对应位置元素相乘，返回矩阵仍然是原型:param arr_w_hat: 3*1 的二维arrayarray([[-1.11069546e-16],[ 8.78503652e-01],[-4.69166570e-02]]):param means_T: 3*1 的二维array:param stds_T: 3*1 的二维array:return: 咱也不知道为啥这么缩放，抄就完事了。。。。标准化的公式： data = (data - data.mean()) / data.std()转化：1. temp = y的均值 * w / y的标准差2. 转化的b = (b - sum(temp)) * y的标准差 + y的均值3. 转化的w = w * y的标准差 / x的标准差最后把w_hat_T恢复成 1 * 3返回"""# data = (data - data.mean()) / data.std()temp = means_T[:-1] * arr_w_hat_T[1:] / stds_T[:-1]arr_w_hat_T[0] = (arr_w_hat_T[0] - np.sum(temp)) * stds_T[-1] + means_T[-1]arr_w_hat_T[1:] = arr_w_hat_T[1:] * stds_T[-1] / stds_T[:-1]return arr_w_hat_T.reshape(1, -1)

注意点：

1、np

`np.matrix(a)`

1、代码里基本运算全都是matrix即二维数据运算，比如w, x全都要一维转化为二维np.matrix()

a = [1,2,3]np.matrix(a)Out[52]: matrix([[1, 2, 3]])a = np.array([1,2,3])np.matrix(a)Out[54]: matrix([[1, 2, 3]])

2、对于dataframe类型的数据，转matrix

data_x_hat = data.iloc[:, 0: cols - 1] # 取dataframe的x：前3列data_y = data.iloc[:, [cols - 1]] # 取dataframe 的y# 获取matrix类型的所有样本X和真实标记y_matrixX = np.matrix(data_x_hat) # 47 * 3y_matrix = np.matrix(data_y) # 47 * 1

矩阵运算

1、星乘：同型矩阵相乘，对应位置的元素相乘，结果仍然同型

（1）直接乘法符号

a = np.array([[1,2,3], [12,3,4]])b = np.array([[1,2,3], [12,3,4]])a * bOut[40]: array([[ 1, 4, 9],[144, 9, 16]])

（2）np.multiply(a, b)

np.multiply(a, b)Out[41]: array([[ 1, 4, 9],[144, 9, 16]])

点乘：内部做内积:`np.dot(a,b)`

A4∗3∗B3∗4→C4∗4A_{4*3} * B_{3*4} \to C_{4*4}A4∗3∗B3∗4→C4∗4

a = np.array([[1,2,3], [12,3,4], [2,3,4], [1,2,3]])b = np.array([[1,2,3,3], [12,3,4,3], [2,3,4,2]])np.dot(a,b)Out[61]: array([[31, 17, 23, 15],[56, 45, 64, 53],[46, 25, 34, 23],[31, 17, 23, 15]])

2、矩阵内元素做n次方：`np.power(a, n)`

np.power(a, 2)Out[45]: array([[ 1, 4, 9],[144, 9, 16]])

构造矩阵和矩阵的形状

1、零矩阵：`np.zeros((m,n), dtype = )`

参数1：形状，元组形式传入参数dtype：元素类型，默认为float

np.zeros((2,3))Out[47]: array([[0., 0., 0.],[0., 0., 0.]])np.zeros((2,3), dtype = int)Out[48]: array([[0, 0, 0],[0, 0, 0]])

2、矩阵转置：`a.T`

aOut[49]: array([[ 1, 2, 3],[12, 3, 4]])a.TOut[50]: array([[ 1, 12],[ 2, 3],[ 3, 4]])

查看矩阵是几行几列：`a.shape`

a = np.array([[1,2,3], [12,3,4]])a.shapeOut[37]: (2, 3)

转换矩阵形状：`a.reshape(m, n)`

如：np.matrixw = [[-1.11069546e-16 8.78503652e-01 -4.69166570e-02]]

a/b = -1表示不关心行数/列数w.reshape(-1,1) 表示把1行3列的matrix变成1列的matrix，形式为

matrix([[-1.11069546e-16],[ 8.78503652e-01],[-4.69166570e-02]])

2、pandas

read_csv获取dataframe

1、除了path不给任何参数：默认head是第一行数据

2、header = 1/2/3…：指定行x为head，数据从head之后开始

data1 = pd.read_csv(path, header = 1)

3、names参数：给每列指定head名

data = pd.read_csv(path, names=["Sizes", "Bedrooms", "Prices"])

4、data.head()：只输出五行（不包括head）

data.xx().values：获取mean,std,min,max等数据的描述值

data.iloc[p1, p2]

对于dataframe对象使用iloc[p1, p2],输出形式都是：dataframe的形式

p1:表示取哪些行：

a:b：提取行a~b-1[c]：提取行c

p2:表示取哪些列

a:b：提取列a~b-1[c]：提取列c

例：

cols = 4data_x_hat = data.iloc[:, 0: cols - 1] # 取dataframe的x：前3列data_y = data.iloc[:, [cols - 1]] # 取dataframe 的y

data.insert()

data.insert(0, 'Ones', 1) # 在第0列，插入一列属性值全为1的列，列名Ones

完整代码

#!/usr/bin/env python# -*- coding: utf-8 -*-# @Time : /4/5 20:52# @Author : cc# @File : multi-variable-gd.py# @Software: PyCharm# 房价预测。 ex1data2.txt：面积、卧室数、房价import numpy as npimport pandas as pdfrom matplotlib import pyplot as pltfrom mpl_toolkits.mplot3d import Axes3D# 代价函数，返回该模型w，b参数下的代价def cost_function(w_hat_matrix, X_matrix, y_matrix) -> float:""" 此方法不管几维都固定不变:param w_hat_matrix: 二维matrix，1 * 3:param X_matrix: 二维matrix，47 * 3:param y_matrix: 二维matrix，47 * 1:return: 代价"""m = len(X_matrix) # 样本数return np.sum(np.power(X_matrix * w_hat_matrix.T - y_matrix, 2)) / (2 * m)# 梯度下降函数def gradientDescent(alpha, iters, w_hat_matrix, X_matrix, y_matrix) -> tuple:""":param alpha: 梯度下降学习率/步长:param iters: 梯度下降次数:param w_hat_matrix: 二维matrix，1 * 3；一般初始设置为[[0, 0, 0]]，初值对梯度下降收敛速度影响大:param X_matrix: 二维matrix，47 * 3。就是公式推导里记的那个超大矩阵X.存放扩展了1的所有样本:param y_matrix: 二维matrix，47 * 1。存放真实标记:return:"""parameters = int(w_hat_matrix.shape[1]) # 3个参数m = len(X_matrix) # 样本数47cur_w_mat_matrix = np.matrix(np.zeros(w_hat_matrix.shape)) # 暂存每次迭代得到的w_hat1*3cost = np.zeros(iters) # 记录每次迭代后的新的代价cost# 迭代iters次for i in range(iters):# 所有样本x_i预测输出和真实标记 y_i的误差。第i行记录了样本x_i的误差error_matrix = X_matrix * w_hat_matrix.T - y_matrix # 97*1# 梯度下降公式：得到新的w_hat：j指向每一列，进行更新（w_1, w_2,..,w_j,...w_d, b）for j in range(parameters):# np.multiply：m*n和m*n的相同下标元素元素相乘，结果还是m*n的矩阵。可以看成多个数乘s = np.sum(np.multiply(error_matrix, X_matrix[:, j])) # 公式：error矩阵和x的第j列相乘cur_w_mat_matrix[0, j] = cur_w_mat_matrix[0, j] - alpha * s / mw_hat_matrix = cur_w_mat_matrix # 每次下降更新w_hatcost[i] = cost_function(w_hat_matrix, X_matrix, y_matrix) # 更新代价，用于观察梯度下降的代价变化曲线return w_hat_matrix, cost# 画结果图def plot_res(transform_w_hat, data_origin, mins, maxs):""":param transform_w_hat: 符合缩放前原始样例的 b,w1,w2,...wm二维array：[[88307.21151185 138.22534685 -7709.05876589]]:param data_origin: 原始数据每列的值：也是二维array形式属性x1列的所有值array：data_origin[:, 0]属性x2... data_origin[:, 1]:param mins: 原始数据每列最小值array: [ 8521 169900]:param maxs: 原始数据每列最大值array: [ 44785 699900]:return: 画出拟合模型图和散点图"""# 建立三维模型fig = plt.figure() # Create a new figure, or activate an existing figure.ax = Axes3D(fig, auto_add_to_figure=False) # 三维坐标轴fig.add_axes(ax) # 给图形fig添加坐标轴 an Axes to the figure.# 设置三维图角度ax.view_init(elev=25, azim=125) # 10 80观察更好# 设置三根轴的名称ax.set_xlabel('Size')ax.set_ylabel('Bedrooms')ax.set_zlabel('Prices')# 设置x1 x2轴范围x1 = np.arange(mins[0], maxs[0] + 1, 1) # x1轴的范围：步长为1x2 = np.arange(mins[1], maxs[1] + 1, 1) # x2轴的范围：步长为1x1, x2 = np.meshgrid(x1, x2) # 生成网格点坐标矩阵，这句话必须有# 画线性回归模型：平面图b, w1, w2 = transform_w_hat[0, 0], transform_w_hat[0, 1], transform_w_hat[0, 2]# 获取系数f = b + w1 * x1 + w2 * x2# 模型：映射关系ax.plot_surface(x1, x2, f, rstride=1, cstride=1, color='red') # 创建平面图（模型）# 创画样例散点图ax.scatter(data_origin[:, 0], data_origin[:, 1], data_origin[:, 2])plt.show()# 每次梯度下降的代价变化图def plot_cost(cost, iters: int):""":param cost: 一维array，第i个元素存放第i次梯度下降时的代价:param iters: 迭代次数，固定1000次:return:"""# 二维坐标轴直接Plt：设置坐标轴名称和标题plt.xlabel("iterations")plt.ylabel("Cost")plt.title("Error vs Traning Epoch")# 画直线，x变化范围为0~迭代次数，y为每次的代价plt.plot(range(iters), cost, color='red')plt.show()# 把数据经过特征缩放（均值标准化）的w_hat变成符合原始数据的w_hat_transformdef w_hat_transform(arr_w_hat_T, means_T, stds_T):""" 同型array相乘相除：对应位置元素相乘，返回矩阵仍然是原型:param arr_w_hat: 3*1 的二维arrayarray([[-1.11069546e-16],[ 8.78503652e-01],[-4.69166570e-02]]):param means_T: 3*1 的二维array:param stds_T: 3*1 的二维array:return: 咱也不知道为啥这么缩放，抄就完事了。。。。标准化的公式： data = (data - data.mean()) / data.std()转化：1. temp = y的均值 * w / y的标准差2. 转化的b = (b - sum(temp)) * y的标准差 + y的均值3. 转化的w = w * y的标准差 / x的标准差最后把w_hat_T恢复成 1 * 3返回"""# data = (data - data.mean()) / data.std()temp = means_T[:-1] * arr_w_hat_T[1:] / stds_T[:-1]arr_w_hat_T[0] = (arr_w_hat_T[0] - np.sum(temp)) * stds_T[-1] + means_T[-1]arr_w_hat_T[1:] = arr_w_hat_T[1:] * stds_T[-1] / stds_T[:-1]return arr_w_hat_T.reshape(1, -1)if __name__ == '__main__':""" pandas得到dataframe类型的数据data SizeBedroomsPrice 3列"""path = 'ex1data2.txt'data = pd.read_csv(path, names=["Sizes", "Bedrooms", "Prices"])# 获取原始数据data的一些描述data_origin = data.values # 二维array，每个子array存放样例的数据means = data.mean().values # 一维array，第i个元素：第i列的均值stds = data.std().values # 一维array，第i个元素：第i列的标准差mins = data.min().values # 一维array，第i个元素：第i列的最小值maxs = data.max().values # 一维array，第i个元素：第i列的最大值'''data特征缩放：均值标准化，让不同特征值差异不要太大，否则梯度下降收敛会很慢'''data = (data - data.mean()) / data.std()'''添加列，用于获得x_hat组成的X：详情见公式推导，为了计算省略bOnes SizeBedroomsPrice 4列'''data.insert(0, 'Ones', 1) # 在第0列，插入一列属性值全为1的列，列名Onescols = data.shape[1] # 列数4"""对于dataframe对象使用iloc[p1, p2]p1:表示取哪些行： a:b：提取行a~b-1[c]：提取行cp2:表示取哪些列a:b：提取列a~b-1[c]：提取列c输出形式都是：dataframe的形式"""data_x_hat = data.iloc[:, 0: cols - 1] # 取dataframe的x：前3列data_y = data.iloc[:, [cols - 1]] # 取dataframe 的y# 获取matrix类型的所有样本X和真实标记y_matrix：有的计算只能二维matrix做X = np.matrix(data_x_hat) # 47 * 3y_matrix = np.matrix(data_y) # 47 * 1w_hat_matrix = np.matrix([0, 0, 0]) # 1*3 初始w一般设置为0# 看看获取的matrix是几行几列： (47, 3) (47, 1) (1, 3)# print(X.shape, y_matrix.shape, w_hat_matrix.shape)# 设置步长和下降次数alpha = 0.01 # 常考率 0.01， 0.03， 0.1， 0.3， 1， 3， 10iters = 1000 # 迭代次数# 梯度下降iters次后，我们获得的res_w_hat能使得代价函数获得【局部最小值】res_w_hat, cost = gradientDescent(alpha, iters, w_hat_matrix, X, y_matrix)print("res_w_hat：", res_w_hat)# 用能使得代价函数获得【局部最小值】的res_w_hat，带入计算局部最小代价res_cost = cost_function(res_w_hat, X, y_matrix)print("res_cost：", res_cost)"""之前获得的res_w_hat，是数据缩放后得到的res_w_hat最终我们要把res_w_hat恢复成与原始数据对应的transform_w_hatdata = (data - data.mean()) / data.std()""""""reshape(a, b)如：np.matrix w = [[-1.11069546e-16 8.78503652e-01 -4.69166570e-02]]a/b = -1表示不关心行数/列数w.reshape(-1,1) 表示把1行3列的matrix变成1列的matrix，形式为matrix([[-1.11069546e-16],[ 8.78503652e-01],[-4.69166570e-02]])pandas变成excel就是只有一列 """# 咋说呢，就是做了个转置res_w_hat_reshape = np.array(res_w_hat.reshape(-1, 1)) # 把1*3的二维matrix变成 3*1 的二维arraymeans_reshape = means.reshape(-1, 1) # 3 * 1 的二维arraystds_reshape = stds.reshape(-1, 1) # 3*1 的二维arraytransform_w_hat = w_hat_transform(res_w_hat_reshape, means_reshape, stds_reshape)print("transform_w_hat:", transform_w_hat)plot_res(transform_w_hat, data_origin, mins, maxs) # 模型plot_cost(cost, iters) # 代价函数曲线