4. 吴恩达机器学习--偏差与方差

利用水库水位变化预测大坝出水量，数据集为：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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

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()
1
2
3
4
5
6

3. 定义带正则化的损失函数

  很明显，图中数据应该是一个非线性模型，但为了展示欠拟合的效果，现假设模型为线性模型：$$h_{\theta }(x)={\theta }_0+{\theta }_{1}x$$
  损失函数为：$$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$$

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))
1
2
3
4
5
6

4. 参数初始化

theta = np.ones(X_train.shape[1])
lamda = 1
reg_cost(theta, X_train, y_train, lamda)
1
2
3

  计算所得的初始损失值为：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)
1
2
3
4
5
6
7
8

  执行梯度下降算法之后，得到的参数 $\theta =[-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
1
2
3
4
5
6
7
8
9
10

  首先介绍一下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)
1
2
3

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()
1
2
3
4
5
6

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
1
2
3
4

2. 计算均值和标准差

def get_means_stds(X):
    means = np.mean(X, axis=0)
    stds = np.std(X, axis=0)
    return means, stds
1
2
3
4

3. 特征归一化处理

def feature_normalize(X, means, stds):
    X[:, 1:] = (X[:, 1:] - means[1:]) / stds[1:]
    
    return X
1
2
3
4

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)
1
2
3
4
5
6
7
8
9
10

5. 训练模型

# lamda过小，则会导致过拟合现象
theta_fit = train_model(X_train_norm, y_train, lamda=0)
1
2

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()
1
2
3
4
5
6
7
8
9
10
11
12

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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

  当训练集误差和验证集误差相近时：偏差/欠拟合
  当验证集误差远大于训练集误差时：方差/过拟合

2. 绘制不同lamda的学习曲线

plot_learning_curve(X_train_norm, y_train, X_val_norm, y_val, lamda=0)
1

plot_learning_curve(X_train_norm, y_train, X_val_norm, y_val, lamda=1)
1

plot_learning_curve(X_train_norm, y_train, X_val_norm, y_val, lamda=100)
1

  当正则项系数较小时，对高次数参数的惩罚力度比较大，容易出现欠拟合现象；当正则项系数较大时，对高次数参数的惩罚力度比较小，容易出现过拟合现象。所以要选取合适的正则项系数。

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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

&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)
1
2
3
4
5

&emps;&emps;此时的测试集误差为：4.397616335103924