利用水库水位变化预测大坝出水量,数据集为:ex5data1.mat。
1. 欠拟合展示
1. 导入库,加载数据集
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from scipy.optimize import minimize
data = loadmat('../data/ex5data1.mat')
# 训练集
X_train, y_train = data['X'], data['y']
# 验证集
X_val, y_val = data['Xval'], data['yval']
# 测试集
X_test, y_test = data['Xtest'], data['ytest']
X_train = np.insert(X_train, 0, 1, axis=1) # 在数据集的第一列插入一列全1
X_val = np.insert(X_val, 0, 1, axis=1)
X_test = np.insert(X_test, 0, 1, axis=1)
2. 对训练集数据进行可视化展示
def plot_data():
fig, ax = plt.subplots()
ax.scatter(X_train[:, 1], y_train)
ax.set(xlabel='change in water level(x)', ylabel='water flowing out og the dam(y)')
plot_data()
3. 定义带正则化的损失函数
很明显,图中数据应该是一个非线性模型,但为了展示欠拟合的效果,现假设模型为线性模型:
h
θ
(
x
)
=
θ
0
+
θ
1
x
h_\theta(x) = \theta_0 + \theta_1x
hθ(x)=θ0+θ1x
损失函数为:
J
(
θ
)
=
1
2
m
∑
i
=
1
m
(
h
θ
(
x
(
i
)
)
−
y
(
i
)
)
2
+
λ
2
m
∑
j
=
1
m
θ
j
2
J(\theta) = \frac{1}{2m}\sum_{i=1}^m{(h_\theta(x^{(i)})-y^{(i)})}^2 + \frac{\lambda}{2m}\sum_{j=1}^m{\theta_j^2}
J(θ)=2m1i=1∑m(hθ(x(i))−y(i))2+2mλj=1∑mθj2
def reg_cost(theta, X, y, lamda):
# 由于这里的y为二维数据,所以需要先将y展开为一维数据
cost = np.sum(np.power(X @ theta - y.flatten(), 2))
reg = theta[1:] @ theta[1:] * lamda
return (cost + reg) / (2 * len(X))
4. 参数初始化
theta = np.ones(X_train.shape[1])
lamda = 1
reg_cost(theta, X_train, y_train, lamda)
计算所得的初始损失值为:303.9931922202643
5. 定义带正则项的梯度下降函数
梯度下降的公式为:
def reg_gradient(theta, X, y, lamda):
grad = (X @ theta - y.flatten()) @ X
reg = lamda * theta
reg[0] = 0
return (grad + reg) / len(X)
reg_gradient(theta, X_train, y_train, lamda)
执行梯度下降算法之后,得到的参数 θ = [ − 15.30301567 , 598.25074417 ] \theta= [-15.30301567, 598.25074417] θ=[−15.30301567,598.25074417]
6. 设计训练过程
def train_model(X, y, lamda):
theta = np.ones(X.shape[1]) # 2*1的全1矩阵
res = minimize(fun=reg_cost,
x0 = theta,
args = (X, y, lamda),
method = 'TNC',
jac = reg_gradient)
return res.x
首先介绍一下minimize()函数:用来优化需要求解的参数。函数的返回值中fun是最优函数值,x是最优自变量
res=opt.minimize(fun, x0, args=(), method=None, jac=None, hess=None,
hessp=None, bounds=None, constraints=(), tol=None,
callback=None, options=None)
fun:该参数就是你要去最小化的损失函数。第一个参数必须为theta且其shape必须为(n,)即一维数组
x0:猜测的初始值
args=():优化的附加参数,默认从第二个开始
method:该参数代表采用的方式,默认是BFGS, L-BFGS-B, SLSQP中的一种,可选TNC
BFGS:拟牛顿法,只能处理无约束优化问题,需要使用一阶导数函数。BFGS 算法性能良好,是无约束优化问题的默认算法。
L-BFGS-B:改进的 BFGS 拟牛顿法,L- 指有限内存,-B 指边界约束,可以处理边界约束条件,需要使用一阶导数函数。L-BFGS_B 算法性能良好,消耗内存量很小,适合处理大规模问题,是边界约束优化问题的默认算法。
SLSQP:序贯最小二乘规划算法,可以处理边界约束、等式约束和不等式约束条件。SLSQP 算法性能良好,是带有约束条件优化问题的默认算法。
TNC:截断牛顿法,可以处理边界约束条件
jac:该参数就是计算梯度的函数,和fun参数类似,第一个必须为theta且其shape必须为(n,)即一维数组,最后返回的梯度也必须为一个一维数组。
7. 计算最佳参数,并对模型可视化
theta_final = train_model(X_train, y_train, lamda=0)
theta_final # [13.08790398, 0.36777923]
plot_data()
plt.plot(X_train[:, 1], X_train @ theta_final, c='r')
plt.show()
2. 过拟合展示
1. 定义特征映射函数
def poly_feature(X, power):
for i in range(2, power+1):
X = np.insert(X, X.shape[1], np.power(X[:, 1], i), axis=1)
return X
2. 计算均值和标准差
def get_means_stds(X):
means = np.mean(X, axis=0)
stds = np.std(X, axis=0)
return means, stds
3. 特征归一化处理
def feature_normalize(X, means, stds):
X[:, 1:] = (X[:, 1:] - means[1:]) / stds[1:]
return X
4. 对数据进行特征映射,并做归一化处理
power = 6
X_train_poly = poly_feature(X_train, power)
X_val_poly = poly_feature(X_val, power)
X_test_poly = poly_feature(X_test, power)
train_means, train_stds = get_means_stds(X_train_poly)
X_train_norm = feature_normalize(X_train_poly, train_means, train_stds)
X_val_norm = feature_normalize(X_val_poly,train_means,train_stds)
X_test_norm = feature_normalize(X_test_poly,train_means,train_stds)
5. 训练模型
# lamda过小,则会导致过拟合现象
theta_fit = train_model(X_train_norm, y_train, lamda=0)
6. 可视化展示
def plot_poly_fit():
plot_data()
x = np.linspace(-60, 60, 100)
xx = x.reshape(100,1)
xx = np.insert(xx, 0, 1, axis=1)
xx = poly_feature(xx, power)
xx = feature_normalize(xx, train_means, train_stds)
plt.plot(x, xx @ theta_fit, 'r--')
plot_poly_fit()
3. 绘制学习曲线
1. 设计绘制学习曲线函数
def plot_learning_curve(X_train, y_train, X_val, y_val, lamda):
x = range(1, len(X_train) + 1)
training_cost = []
cv_cost = []
for i in x:
res = train_model(X_train[:i, :], y_train[:i, :], lamda)
training_cost_i = reg_cost(res, X_train[:i, :], y_train[:i, :], lamda)
cv_cost_i = reg_cost(res, X_val, y_val, lamda)
training_cost.append(training_cost_i)
cv_cost.append(cv_cost_i)
plt.plot(x, training_cost, label='training cost')
plt.plot(x, cv_cost, label='cv cost')
plt.legend()
plt.xlabel('number of training examples')
plt.ylabel('error')
plt.show()
当训练集误差和验证集误差相近时:偏差/欠拟合
当验证集误差远大于训练集误差时:方差/过拟合
2. 绘制不同lamda的学习曲线
plot_learning_curve(X_train_norm, y_train, X_val_norm, y_val, lamda=0)
plot_learning_curve(X_train_norm, y_train, X_val_norm, y_val, lamda=1)
plot_learning_curve(X_train_norm, y_train, X_val_norm, y_val, lamda=100)
当正则项系数较小时,对高次数参数的惩罚力度比较大,容易出现欠拟合现象;当正则项系数较大时,对高次数参数的惩罚力度比较小,容易出现过拟合现象。所以要选取合适的正则项系数。
3. 绘制误差函数随lamda变化的曲线
lamdas = [0,0.001,0.003,0.01,0.03,0.1,0.3,1,3,10]
training_cost = []
cv_cost = []
for lamda in lamdas:
res = train_model(X_train_norm,y_train,lamda)
tc = reg_cost(res,X_train_norm,y_train,lamda = 0)
cv = reg_cost(res,X_val_norm,y_val,lamda=0)
training_cost.append(tc)
cv_cost.append(cv)
plt.plot(lamdas,training_cost,label='training cost')
plt.plot(lamdas,cv_cost,label='cv cost')
plt.legend()
plt.xlabel('lamda')
plt.ylabel('cost')
plt.show()
&emps;&emps;由图像可找出使得验证集误差最小的
λ
\lambda
λ
lamdas[np.argmin(cv_cost)] # lamda = 3
res = train_model(X_train_norm,y_train,lamda =3)
test_cost = reg_cost(res,X_test_norm,y_test,lamda = 0)
print(test_cost)
&emps;&emps;此时的测试集误差为:4.397616335103924