一：梯度下降法推导公式：
二：当训练集为1维时
#进行数据分析所需库，可以看做是对numpy工具的补充 import pandas as pd import numpy as np #应该把Seaborn视为matplotlib的补充,作图所用工具，在大多数情况下使用seaborn就能做出很具有吸引力的图，而使用matplotlib就能制作具有更多特色的图 import seaborn as sns import matplotlib.pyplot as plt #设置绘画的图标格式和颜色 sns.set(context="notebook", style="whitegrid", palette="dark") #读取数据并赋予列名，此时df共有两列，列名分别为population和profit df = pd.read_csv('ex1data1.txt', names=['population', 'profit']) def get_x(df): #ones:格式化为一个m列的全为1的向量 ones = pd.DataFrame({'ones':np.ones(len(df))}) #将date格式化为ones向量右边加上df矩阵的矩阵 data = pd.concat([ones, df], axis=1) #返回data矩阵的前两列，返回m * 2的矩阵 return data.iloc[:, :-1].as_matrix() # 这个操作返回 ndarray,不是矩阵 def get_y(df): #返回df矩阵的最后一列，返回m × 1的向量 return np.array(df.iloc[:, -1]) def normalize_feature(df): # 特征缩放,对df中的两列数据分别进行特征缩放 ：（x - x平均值）/ x方差 return df.apply(lambda column: (column - column.mean())/column.std()) X = get_x(df) print(X.shape, type(X)) #X为m × 2的矩阵 y = get_y(df) print(y.shape, type(y)) #y为m × 1的向量 theta = np.zeros(X.shape[1])#X.shape[1]=2,代表特征数n X.shape(0) = 97;X.shape(1) = 2;将theta初始化为2行的零向量（列向量） def lr_cost(theta, X, y): # """ # 计算theta固定时，此时的代价函数值 # X: R(m*n), m 样本数, n 特征数 # y: R(m) # theta : R(n), 线性回归的参数 # """ m = X.shape[0]#m为样本数 #X矩阵和theta矩阵的点积所得矩阵（m × 1） - 矩阵y（m × 1） inner = X @ theta - y # R(m*1)，X @ theta等价于X.dot(theta) #square_sum为inner的每个元素平方之和，即二范数 square_sum = inner.T @ inner #得到此时代价函数的值 cost = square_sum / (2 * m) return cost #将代价函数对theta求导，即求梯度 def gradient(theta, X, y): #m为样本数 m = X.shape[0] #这里的inner为代价函数对theta求导的结果 inner = X.T @ (X @ theta - y) # (m,n).T @ (m, 1) -> (n, 1)，X @ theta等价于X.dot(theta) return inner / m #梯度下降函数 def batch_gradient_decent(theta, X, y, epoch, alpha=0.01): # 拟合线性回归，返回参数和代价 # epoch: 批处理的轮数 # alpha: theta移动的步长 # """ #得到初始theta时的代价 cost_data = [lr_cost(theta, X, y)] # 拷贝一份，不和原来的theta混淆 _theta = theta.copy() #开始迭代epoch次 for _ in range(epoch): #根据当前梯度更新theta _theta = _theta - alpha * gradient(_theta, X, y) #记录此时的代价 cost_data.append(lr_cost(_theta, X, y)) return _theta, cost_data #迭代500次求最小代价和所对应的theta epoch = 500 final_theta, cost_data = batch_gradient_decent(theta, X, y, epoch) # 计算最终的代价 lr_cost(final_theta, X, y) #画出代价函数值变化图 #可以看到从第二轮代价数据变换很大，接下来平稳了 ax = sns.tsplot(cost_data, time=np.arange(epoch+1)) ax.set_xlabel('epoch') ax.set_ylabel('cost') plt.show() b = final_theta[0] # intercept，Y轴上的截距 m = final_theta[1] # slope，斜率 #画出原数据点和线性回归的最终结果 plt.scatter(df.population, df.profit, label="Training data") plt.plot(df.population, df.population*m + b, label="Prediction") plt.legend(loc=2) plt.show()
三：当训练集为多维时
当训练集为多维时只需将上述代码稍作改动，完整代码如下：
#进行数据分析所需库，可以看做是对numpy工具的补充 import pandas as pd import numpy as np #应该把Seaborn视为matplotlib的补充,作图所用工具，在大多数情况下使用seaborn就能做出很具有吸引力的图，而使用matplotlib就能制作具有更多特色的图 import seaborn as sns import matplotlib.pyplot as plt #设置绘画的图标格式和颜色 sns.set(context="notebook", style="whitegrid", palette="dark") #读取数据并赋予列名，此时df共有两列，列名分别为population和profit df = pd.read_csv('ex1data2.txt', names=['square', 'bedrooms', 'price']) df.head() def get_x(df): #ones:格式化为一个m列的全为1的向量 ones = pd.DataFrame({'ones':np.ones(len(df))}) #将date格式化为ones向量右边加上df矩阵的矩阵 data = pd.concat([ones, df], axis=1) #返回data矩阵的前两列，返回m * 2的矩阵 return data.iloc[:, :-1].as_matrix() # 这个操作返回 ndarray,不是矩阵 def get_y(df): #返回df矩阵的最后一列，返回m × 1的向量 return np.array(df.iloc[:, -1]) def normalize_feature(df): # 特征缩放,对df中的两列数据分别进行特征缩放 ：（x - x平均值）/ x方差 return df.apply(lambda column: (column - column.mean())/column.std()) data = normalize_feature(df) #特征缩放 X = get_x(data) print(X.shape, type(X)) #X为m × 2的矩阵 y = get_y(data) print(y.shape, type(y)) #y为m × 1的向量 theta = np.zeros(X.shape[1])#X.shape[1]=2,代表特征数n X.shape(0) = 97;X.shape(1) = 2;将theta初始化为2行的零向量（列向量） def lr_cost(theta, X, y): # # 计算theta固定时，此时的代价函数值 # X: R(m*n), m 样本数, n 特征数 # y: R(m) # theta : R(n), 线性回归的参数 # m = X.shape[0]#m为样本数 #X矩阵和theta矩阵的点积所得矩阵（m × 1） - 矩阵y（m × 1） inner = X @ theta - y # R(m*1)，X @ theta等价于X.dot(theta) #square_sum为inner的每个元素平方之和，即二范数 square_sum = inner.T @ inner #得到此时代价函数的值 cost = square_sum / (2 * m) return cost #将代价函数对theta求导，即求梯度 def gradient(theta, X, y): #m为样本数 m = X.shape[0] #这里的inner为代价函数对theta求导的结果 inner = X.T @ (X @ theta - y) # (m,n).T @ (m, 1) -> (n, 1)，X @ theta等价于X.dot(theta) return inner / m #梯度下降函数 def batch_gradient_decent(theta, X, y, epoch, alpha=0.01): # 拟合线性回归，返回参数和代价 # epoch: 批处理的轮数 # alpha: theta移动的步长 # #得到初始theta时的代价 cost_data = [lr_cost(theta, X, y)] # 拷贝一份，不和原来的theta混淆 _theta = theta.copy() #开始迭代epoch次 for _ in range(epoch): #根据当前梯度更新theta _theta = _theta - alpha * gradient(_theta, X, y) #记录此时的代价 cost_data.append(lr_cost(_theta, X, y)) return _theta, cost_data base = np.logspace(-1, -5, num=4) candidate = np.sort(np.concatenate((base, base*3))) epoch = 50 fig, ax = plt.subplots(figsize=(16, 9)) for alpha in candidate: _, cost_data = batch_gradient_decent(theta, X, y, epoch, alpha=alpha) ax.plot(np.arange(epoch+1), cost_data, label=alpha) ax.set_xlabel('epoch', fontsize=18) ax.set_ylabel('cost', fontsize=18) ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.) ax.set_title('learning rate', fontsize=18) plt.show() ''' #画出代价函数值变化图 #可以看到从第二轮代价数据变换很大，接下来平稳了 ax = sns.tsplot(cost_data, time=np.arange(epoch+1)) ax.set_xlabel('epoch') ax.set_ylabel('cost') plt.show() '''