赞
踩
在练习的前半部分,您将实现正则化线性回归,利用水库水位的变化来预测从大坝流出的水量。在下半部分中,您将对调试学习算法进行一些诊断,并检查偏差和偏差的影响。
首先,我们将可视化包含水位变化的历史记录的数据集,X和从大坝流出的水量y。这个数据集包括三个部分
训练集X,y,用于训练模型
验证集Xval,yval,用于选择正则化参数
测试集Xtest,ytest, 用于评估性能
首先我们需要先导入所用到的库
- import numpy as np
- import scipy.io as sio
- import scipy.optimize as opt
- import pandas as pd
- import matplotlib.pyplot as plt
- import seaborn as sns
接着导入本次实验中所用到的数据,在这里我们可以将数据集分别提出来,并查看其大小
- # 读入数据
- path = r'E:\Code\ML\ml_learning\ex5-bias vs variance\ex5data1.mat'
- data = sio.loadmat(path)
- X, y, Xval, yval, Xtest, ytest = data['X'], data['y'], data['Xval'], data['yval'], data['Xtest'], data['ytest']
- # X.shape, y.shape, Xval.shape, yval.shape, Xtest.shape, ytest.shape
- # ((12, 1), (12, 1), (21, 1), (21, 1), (21, 1), (21, 1))
接下来进行数据可视化
- # 可视化
- def plotdata(X,y):
- plt.figure(figsize=(12,8))
- plt.scatter(X, y,c='r', marker='x')
- plt.xlabel('Change in water level (x)')
- plt.ylabel('Water flowing out of the dam (y)')
- plt.grid(True)
- plotdata(X,y)
正则化线性回归的代价函数如下
其中λ是控制正则化程度的正则化参数,因此可以用于防止过拟合,正则化项对总代价J(θ)施加惩罚,随着模型参数θj的增大,惩罚也增大,需要注意的是不惩罚θ0项,即索引从1开始
- # 插入x0 = 1
- X, Xval, Xtest = [np.insert(x, 0, np.ones(x.shape[0]), axis=1) for x in (X, Xval, Xtest)]
- # X.shape, Xval.shape, Xtest.shape
- # ((12, 2), (21, 2), (21, 2))
-
- # 代价函数
- def costReg(theta, X, y, l=1):
- m = X.shape[0] # 12
- theta = np.matrix(theta) # (1,2)
- X = np.matrix(X) # (12,2)
- y = np.matrix(y) # (12,1)
- inner = (X * theta.T) - y # (12,1)
- part1 = float((1 / (2 * m)) * inner.T * inner) #(1,1)
- part2 = float((l / (2 * m)) * theta[:,1] * theta[:,1].T) # (1,1)
- cost = part1 + part2
-
- return cost
其中当初始θ = [1;1]时,所得输出为303.993,在这里一定要注意公式后半部分正则化处不包括θ0这一项,不然输出就变成304.04,反复检查了好几遍才发现这里出了问题。
相应地,正则化线性回归梯度公式定义为
- # 正则化梯度下降
- def gradientReg(theta, X,y, l=1):
- m = X.shape[0] # 12
- theta = np.matrix(theta) # (1,2)
- X = np.matrix(X) # (12,2)
- y = np.matrix(y) # (12,1)
- # (2,12) * ((12,2) * (2,1)) - (12, 1)) = (2,1)
- inner = (1 / m) * X.T * ((X * theta.T) - y)
- reg = (l /m) * theta # (1,2)
- reg[0,0] = 0 #第0项不正则化
- return inner + reg.T # (2,1)
代码中需要注意的就是计算θ0的梯度不需要正则化
当λ = 1所算得初始梯度为[-15.303; 598.250]
在这部分中,我们将λ设为0, 因为目前模型参数较少,正则化不会起什么大的帮助,并调用scipy.optimize中的minimize函数计算参数最优解,并将拟合函数与原始数据画在一起
- # 拟合线性回归
- final_theta = opt.minimize(fun=costReg, x0=theta, args=(X, y, 0), method='TNC', jac=gradientReg, options={'disp': True}).x
- def plotdata1(theta, X, y):
- fig,ax = plt.subplots(figsize=(12,8))
- plt.scatter(X[:,1], y, c='r', label="Training data")
- plt.plot(X[:,1], X @ theta, c = 'b', label="Prediction")
- ax.set_xlabel("water_level")
- ax.set_ylabel("flow")
- ax.legend()
- plt.show()
-
- plotdata1(final_theta, X, y)
最佳拟合直线告诉我们,由于数据具有非线性,模型与数据的拟合不是很好。虽然如图所示可视化最佳拟合是调试学习算法的一种可能方法,但可视化数据和模型并不总是那么容易。于是我们可以实现一个生成学习曲线的函数,即使数据不容易可视化,它也可以帮助调试学习算法。
机器学习中的一个重要概念是偏差-方差权衡。偏差代表着模型与数据的拟合程度,高偏差的模型容易出现拟合不足,而方差代表着模型的泛化能力,高方差的模型对训练数据会出现过度拟合,即对新的样本泛化能力很差。我们可以通过绘制训练集与验证集的误差,来判断偏差与方差问题
为了绘制学习曲线,我们需要不同训练集大小的训练和交叉验证集误差。为了获得不同的训练集大小,你应该使用原始训练集X的不同子集。具体来说,对于一个训练集大小为i的集合,你应该使用前i个示例(即X(1:i,:)和y(1:i))。
模型的最优参数可以通过前面提到的minimize函数进行求解,得出模型参数后开始计算训练集与验证机上的误差,数据集的训练错误定义为
特别需要注意的是,训练误差不包括正则化项。计算训练误差的一种方法是使用您现有的代价函数并设置λ为0,只有当使用它来计算训练误差和交叉验证误差时。当计算训练集误差时,确保你是在训练子集(即X(1:n,:)和y(1:n))上计算它(而不是整个训练集)。但是,对于交叉验证错误,需要在整个交叉验证集上计算它。最后将计算的误差存储在向量、误差序列和误差值中,以便于可视化
- # 学习曲线
- def linear_regression(X,y,l=1):
- """求出最优参数"""
- theta = np.ones(X.shape[1]) # 初始化参数
- # 训练参数
- res = opt.minimize(fun=costReg,x0=theta,args=(X,y,l),method='TNC',jac=gradientReg, options={'disp': True})
- return res.x
-
- def plot_learning_curve(X, y, Xval, yval, l):
- """画出学习曲线"""
- m = X.shape[0]
- training_cost, cv_cost = [], []
-
- for i in range(1,m+1):
- res = linear_regression(X[:i, :], y[:i], l)
-
- tc = costReg(res, X[:i, :], y[:i], 0)
- cv = costReg(res, Xval, yval, 0)
-
- training_cost.append(tc)
- cv_cost.append(cv)
-
- plt.figure(figsize=(12,8))
- plt.plot(np.arange(1, m+1), training_cost, label='training cost')
- plt.plot(np.arange(1, m+1), cv_cost, label='cv cost')
- plt.legend()
- plt.xlabel("Number of training examples")
- plt.ylabel("Error")
- plt.title("Learing curve for linear regression")
- plt.grid(True)
- plt.show()
-
- plot_learning_curve(X, y, Xval, yval, 0)
图中可知随着训练样本增加,两者误差都很高,说明线性回归不能很好拟合数据集,出现了高偏差问题。
在前面的模型中出现的问题是,它对数据过于简单,导致拟合很差,所以采用添加更多特征来解决此问题。对于多项式回归,假设形式如下,增加的特征为原始值的各种幂次。
- def poly_features(X, power):
- """
- 多项式特征
- 每次在X最后一列添加次方项
- 从第二列开始插入,因为X本身含有一列x0 = 1
- """
- 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
增加了特征之后,由于增加的特征为各次幂,故数据差异很大,比如x = 40,它的八次方就达到了10的12次方,所以我们需要对特征向量进行特征标准化。
- def get_means_std(X):
- """获得训练集的均值和误差"""
- means = np.mean(X,axis = 0) # 按列
- # ddof = 1 求样本标准差
- # ddof = 1 求总体标准差
- stds = np.std(X, axis=0, ddof=1)
-
- return means, stds
-
-
- def featureNormalize(X, means, stds):
- """标准化"""
- X_norm = X.copy()
- X_norm[:,1:] = (X_norm[:,1:] - means[1:]) / stds[1:]
- return X_norm
在画学习曲线图之前我们需要进行数据预处理
- # 数据处理
- power = 6 # 扩展到x的6次方
- # 均值与标准差
- train_means, train_stds = get_means_std(poly_features(X, power))
- # 标准化
- X_norm = featureNormalize(poly_features(X, power), train_means, train_stds)
- Xval_norm = featureNormalize(poly_features(Xval, power), train_means, train_stds)
- Xtest_norm = featureNormalize(poly_features(Xtest, power), train_means, train_stds)
接着画出学习曲线,这里画的是当λ = 0的情况
- def plot_fit(means, stds, l):
- """画出拟合曲线"""
- theta = linear_regression(X_norm,y, l)
-
- x = np.linspace(-75,55,50)
-
- xmat = x.reshape(-1, 1) # (50,)->(50,1)
- xmat = np.insert(xmat,0,1,axis=1) # 添加x0 = 1
-
- Xmat = poly_features(xmat, power) # 增加特征
- Xmat_norm = featureNormalize(Xmat, means, stds) # 特征规范化
-
-
- plotdata(X[:,1], y) # 画出原始数据
- plt.plot(x, Xmat_norm @ theta, 'b--')# 画出拟合曲线
-
- plot_fit(train_means, train_stds, 0)
- plot_learning_curve(X_norm, y, Xval_norm, yval, 0) # 画出学习曲线
通过拟合图我们可以看得出多项式拟合能够很好跟随数据。通过学习曲线图可知,随着训练样本增加,训练误差比较低,虽然验证误差整体趋势是下降,但是最后的误差还是很大,训练误差与验证误差存在差距,说明模型对训练数据存在过拟合,高方差的问题,不能很好地推广。
而当λ = 1时,可以看出训练误差和验证误差都收敛到一个比较小的值,代表拟合程度还不错。
而当λ = 100 时,可以看出拟合效果很差,训练误差和验证误差都收敛到一个很大的值,说明正则化过多,模型无法拟合数据。
λ的取值可以选为[0., 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1., 3., 10.]
- lambdas = [0., 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1., 3., 10.]
- errors_train, errors_val = [], []
- for l in lambdas:
- theta = linear_regression(X_norm, y, l)
- errors_train.append(costReg(theta,X_norm,y,0)) # 记得把lambda = 0
- errors_val.append(costReg(theta,Xval_norm,yval,0))
-
- plt.figure(figsize=(8,5))
- plt.plot(lambdas,errors_train,label='Train')
- plt.plot(lambdas,errors_val,label='Cross Validation')
- plt.legend()
- plt.xlabel('lambda')
- plt.ylabel('Error')
- plt.grid(True)
lambdas[np.argmin(errors_val)] # 3.0
可以得出验证误差最小的时候λ = 3
- theta = linear_regression(X_norm, y, 3)
- print('test cost(l={}) = {}'.format(3, costReg(theta, Xtest_norm, ytest, 0))
- # test cost(l=3) = 4.407884454040075
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。