'''在本练习中，您将实现正则化的线性回归和多项式回归，并使用它来研究具有不同偏差-方差属性的模型.在前半部分的练习中，你将实现正则化线性回归，以预测水库中的水位变化，从而预测大坝流出的水量。在下半部分中，您将通过一些调试学习算法的诊断，并检查偏差 v.s. 方差的影响。我们将从可视化数据集开始，其中包含水位变化的历史记录，x，以及从大坝流出的水量，y。这个数据集分为了三个部分：training set 训练集：训练模型cross validation set 交叉验证集：选择正则化参数test set 测试集：评估性能，模型训练中不曾用过的样本'''import numpy as npimport matplotlib.pyplot as pltfrom scipy.io import loadmatimport scipy.optimize as opt# 读取数据path = 'data/ex5data1.mat' # 数据路径data = loadmat(path) # 加载数据X, y = data['X'], data['y'] # 获取 训练集Xval, yval = data['Xval'], data['yval'] # 获取 交叉验证集Xtest, ytest = data['Xtest'], data['ytest'] # 训练集X = np.insert(X, 0, 1, axis=1) # X的第一列插入1Xval = np.insert(Xval, 0, 1, axis=1) # Xval的第一列插入1Xtest = np.insert(Xtest, 0, 1, axis=1) # Xtest的第一列插入1# 打印X，Xval，Xtest的大小# print('X={},y={}'.format(X.shape, y.shape))# print('Xval={},yval={}'.format(Xval.shape, yval.shape))# print('Xtest={},ytest={}'.format(Xtest.shape, ytest.shape))'''结果：X=(12, 2),y=(12, 1)Xval=(21, 2),yval=(21, 1)Xtest=(21, 2),ytest=(21, 1)'''def plotData(): # 将训练集的数据可视化为散点图plt.figure(figsize=(8, 5)) # 图大小plt.scatter(X[:, 1:], y, c='r', marker='x') # X轴，y轴，颜色，点样式plt.xlabel('Change in water level (x)')plt.ylabel('Water flowing out of the dam (y)')plt.grid(True) # 生成网格# plotData() #运行画图函数# plt.show()# 计算带正则化项的线性代价值def costReg(theta, X, y, l):"""不正则化 theta0theta:(n+1,)X:(m, n+1)y:(m, 1)"""cost = ((X @ theta - y.flatten()) ** 2).sum()regterm = l * (theta[1:] @ theta[1:]) # 正则化项return (cost + regterm) / (2 * len(X)) # 返回带正则化项的线性代价值# 使用theta=[1，1]，lambda = 1运行代价计算# theta = np.ones(X.shape[1]) # [1,1]# print(costReg(theta, X, y, l=1)) # 结果：303.9931922202643# 计算带正则化项的线性梯度值def gradientReg(theta, X, y, l):"""不正则化 theta0theta:(n+1,) theta:(2,)X:(m, n+1) X: (12, 2)y:(m, 1) y: (12, 1)返回结果：(2,)"""grad = (X @ theta - y.flatten()) @ Xregterm = l * theta # 正则化项regterm[0] = 0return (grad + regterm) / len(X) # 返回带正则化项的线性梯度值# 使用theta=[1，1]，lambda = 1运行梯度计算# theta = np.ones(X.shape[1])# print(gradientReg(theta, X, y, 1)) # 结果：[-15.30301567 598.25074417]# 拟合线性回归def trainLinearReg(X, y, l):theta = np.ones(X.shape[1]) # 初始化theta为[1，1]res = opt.minimize(fun=costReg, x0=theta, args=(X, y, l), method='TNC', jac=gradientReg)return res.x # 结果返回优化后的theta:(2,)fit_theta = trainLinearReg(X, y, 0)print('优化后的theta:', fit_theta) # 结果：[13.08790367 0.36777923]plotData()plt.plot(X[:, 1:], X @ fit_theta) # 画线性拟合直线'''这里我们把λ = 0，因为线性回归只有两个参数，这么低的维度，正则化并没有用。从图中可以看到，拟合最好的这条直线告诉我们这个模型并不适合这个数据。接下来，您将实现一个函数来生成学习曲线，它可以帮助您调试学习算法。''''''机器学习中一个重要的概念是偏差（bias）和方差（variance）的权衡。高偏差意味着欠拟合，高方差意味着过拟合。在这部分练习中，您将在学习曲线上绘制训练误差和验证误差，以诊断bias-variance问题。''''''训练样本数量从1开始逐渐增加，训练出不同的参数向量θ。接着通过交叉验证样本Xval计算验证误差。切记此时不要使用正则化，将λ=0。计算交叉验证代价时记得整个交叉验证集来计算。'''# 画出学习曲线，即交叉验证误差和训练误差随样本数量的变化的变化def plot_learning_curve(X, y, Xval, yval, l):xx = range(1, len(X) + 1) # [1,2,3,4,...,12]training_cost, cv_cost = [], []for i in xx:res = trainLinearReg(X[:i], y[:i], l) # 样本数量逐层增加,res存储每次的theta值training_cost_i = costReg(res, X[:i], y[:i], 0) # 训练集代价cv_cost_i = costReg(res, Xval, yval, 0) # 交叉验证集代价training_cost.append(training_cost_i)cv_cost.append(cv_cost_i)# 画学习曲线图plt.figure(figsize=(8, 5))plt.plot(xx, training_cost, label='training cost')plt.plot(xx, cv_cost, label='cv cost')plt.legend()plt.xlabel('Number of training examples')plt.ylabel('Cost')plt.title('Learning curve for linear regression')plt.grid(True)plot_learning_curve(X, y, Xval, yval, l=0)# 从图中看出来，随着样本数量的增加，训练误差和交叉验证误差都很高，这属于高偏差，欠拟合。'''我们的线性模型对于数据来说太简单了，导致了欠拟合(高偏差)。接下来，您将通过添加更多的特性来解决这个问题。数据预处理:X，Xval，Xtest都需要添加多项式特征，这里我们选择增加到6次方(不要忘了标准化)。'''def genPolyFeatures(X, power):""" 添加多项式特征:每次在最后一列插入第二列的 2,3,...,power 次方（第一列为偏置）从二次方开始插入（因为本身含有一列一次方）"""Xpoly = X.copy() # 浅拷贝for i in range(2, power + 1):Xpoly = np.insert(Xpoly, Xpoly.shape[1], np.power(Xpoly[:, 1], i), axis=1) # 按列插入return Xpoly # 返回添加多项式特征后的数据集def get_means_std(X):"""获取训练集的均值和误差，用来标准化所有数据。"""means = np.mean(X, axis=0) # 对各列求均值，返回 1* n 矩阵"""numpy.mean(a, axis, dtype, out，keepdims )函数功能：求取均值经常操作的参数为axis，以m * n矩阵举例：axis 不设置值，对 m*n 个数求均值，返回一个实数axis = 0：压缩行，对各列求均值，返回 1* n 矩阵axis =1 ：压缩列，对各行求均值，返回 m *1 矩阵"""stds = np.std(X, axis=0, ddof=1) # 计算每一列的无偏样本标准差"""numpy.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=<no value>)函数功能：求标准差axis=0时，表示求每一列标准差axis=1时，表示求每一行标准差axis=None时，表示求全局标准差。ddof=0时，计算有偏样本标准差(一般在拥有所有数据的情况下，计算所有数据的标准差时使用，即最终除以n)ddof=1时，表示计算无偏样本标准差(最终除以n-1)"""return means, stdsdef featureNormalize(myX, means, stds):"""标准化所有数据集应该都用训练集的均值和样本标准差处理。"""X_norm = myX.copy()#浅拷贝X_norm[:, 1:] = X_norm[:, 1:] - means[1:]X_norm[:, 1:] = X_norm[:, 1:] / stds[1:]return X_norm# 获取添加多项式特征以及标准化之后的数据集。power = 6 # 扩展到x的6次方train_means, train_stds = get_means_std(genPolyFeatures(X, power)) # 获得训练集的均值和标准差X_norm = featureNormalize(genPolyFeatures(X, power), train_means, train_stds) # 扩展特征后的训练集Xval_norm = featureNormalize(genPolyFeatures(Xval, power), train_means, train_stds) # 扩展特征后的交叉验证集Xtest_norm = featureNormalize(genPolyFeatures(Xtest, power), train_means, train_stds) # 扩展特征后的测试集# 画出多项式回归的拟合曲线def plot_fit(means, stds, l):theta = trainLinearReg(X_norm, y, l) # 得到扩大特征后的训练集的最优theta值 (7, )print('此时的theta:', theta)x = np.linspace(-75, 55, 50) # 50表示要生成的样本数xmat = x.reshape(-1, 1) # 重塑成一列 (50,1)xmat = np.insert(xmat, 0, 1, axis=1) # 在第一列插入1 (50,2)Xmat = genPolyFeatures(xmat, power) # 扩大特征 (50,7)Xmat_norm = featureNormalize(Xmat, means, stds) # 标准化plotData() # 画原数据散点图plt.plot(x, Xmat_norm @ theta, 'b--') # 画拟合曲线plot_fit(train_means, train_stds, 0)# 此时的theta: [ 11.21817067 11.36845531 13.44911626 10.73084633 -4.44505759 -11.89800257 -5.06836084]plot_learning_curve(X_norm, y, Xval_norm, yval, 0)plot_fit(train_means, train_stds, 100)# 此时的theta: [11.21758894 0.98876228 0.30025154 0.77917857 0.11521526 0.5903675 -0.0118514 ]plot_learning_curve(X_norm, y, Xval_norm, yval, 100)'''l=0时，不进行正则化，图像显示过拟合特征。l=100时，惩罚过多，图像显示欠拟合特征。'''# 使用交叉验证集来选择最合适的正则化lambda值lambdas = [0., 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1., 3., 10.] # lambda候选值error_train, error_val = [], []for l in lambdas:theta = trainLinearReg(X_norm, y, l)error_train.append(costReg(theta, X_norm, y, 0)) # 这里计算error时，令l=0error_val.append(costReg(theta, Xval_norm, yval, 0))'''正则化是用来进行梯度下降的，惩罚某些特征。当梯度下降完成，用得到的最优参数计算损失函数，就不用正则化那一项了。因此这里计算error时，令l=0。'''l = lambdas[np.argmin(error_val)] # np.argmin(error_val)：获得error_val最小值处的索引号print('最优lambda值为:', l)theta = trainLinearReg(X_norm, y, l)print('此时theta值为：', theta)print('test cost=', costReg(theta, Xtest_norm, ytest, 0))'''运行结果：最优lambda值为: 3.0此时theta值为： [11.21758982 6.77857869 3.98943697 3.99183811 2.24949203 2.45682122 1.25070152]test cost= 4.755261169964735若power=8，则结果为：最优lambda值为: 3.0此时theta值为： [11.2175881 6.65682303 3.89421607 3.75899131 2.25315894 2.20639341 1.33198932 1.35014165 0.76188972]test cost= 3.859893821669393'''# 画lambda-cost图plt.figure(figsize=(8, 5))plt.plot(lambdas, error_train, label='Train')plt.plot(lambdas, error_val, label='Cross Validation')plt.legend()plt.xlabel('lambda')plt.ylabel('Error')plt.grid(True)plt.show()
运行结果：
线性回归及其学习曲线：
多项式回归及其学习曲线：
过拟合：
欠拟合：
lambda的选取：
正则化代价函数：
正则化线性梯度值:
训练集代价函数:
多项式回归:
linspace函数:
reshape函数:
在创建DataFrame的时候常常使用reshape来更改数据的列数和行数。
在使用了reshape（-1，1）之后，数据集变成了一列。reshape(1,-1)变成了一行。
