【图解例说机器学习】模型选择：偏差与方差 (Bias vs. Variance)

机器学习的过程大致分为三步：1）模型假设，比如我们假设模型是线性回归，还是多项式回归，以及其阶数的选择；2）误差函数定义，比如我们假设误差函数是均方误差，还是交叉熵；3）参数求解，比如使用正规方程，还是梯度下降等。

这篇文章主要讨论模型的选择问题，下面以多项式回归为例进行说明

---

一个例子：多项式回归中的阶数选择

在前面的文章【图解例说机器学习】线性回归中，我们定义了广义的线性回归模型，其表达式为：

$$\hat y=\omega_0+\sum\limits_{j=1}^{M}\omega_j\phi_j(\mathrm x)=\omega_0+\mathrm w^{\mathrm T}\phi(\mathrm x)\tag{1}$$

当$D=1,\phi_j(\mathrm x)=x^j$时，公式(1)可以表示为：

$$\hat y=\omega_0+\omega_1x+\omega_2x^2+\cdots+\omega_Mx^M\tag{2}$$

此时，线性回归就变成了$M$阶多项式回归。

当$M$及误差函数给定时，我们就可以通过梯度下降法求解得到$\mathrm w$。但是，$M$的选择对预测的结果影响较大。从公式可以看出$M$越大，模型越复杂，其函数表达式集合包含了$M$取值较小的情况。从这种角度来看，$M$取值越大越好。但是，一般来说训练数据较少，当$M$取值较大时，复杂的模型会过度学习训练数据间的关系，导致其泛化能力较差。

---

这里我们通过一个实例来形象化$M$对算法的影响。这里我们假设实际的函数表达式为

$$y=\sin(2\pi x)+\epsilon\tag{3}$$

其中，$\epsilon$是一个高斯误差值。通过公式(3)我们产生10个样例点$(x_i,y_i)$。在给定不同$M$值时，我们使用正规方程法或梯度下降法可以得到最佳的函数表达式。在这里，我们采用正规方程法 (见【图解例说机器学习】线性回归中公式(12))，得到最优参数：

$$\mathrm{\bar w}=[\bar\phi^{\mathrm T}(\mathrm X)\bar\phi(\mathrm X)]^{-1}\bar\phi^{\mathrm T}(\mathrm X)\mathrm y\tag{4}$$

其中，这里的$\bar{\phi}^{\mathrm T}(\mathrm X)$根据公式(2)和【图解例说机器学习】线性回归中的公式(12)可得

$$\bar\phi(\mathrm X)= \left\{\begin{matrix} \phi_0(\mathrm x_1) & \phi_1(\mathrm x_1) & \cdots & \phi_M(\mathrm x_1)\\ \phi_0(\mathrm x_2) & \phi_1(\mathrm x_2) & \cdots & \phi_M(\mathrm x_2)\\ \vdots & \vdots & \cdots &\vdots \\ \phi_0(\mathrm x_N) & \phi_1(\mathrm x_N) & \cdots & \phi_M(\mathrm x_N) \end{matrix} \right\}= \left\{\begin{matrix} 1 & \mathrm x_1^1 & \cdots & \mathrm x_1^{M}\\ 1 & \mathrm x_2^1 & \cdots & \mathrm x_2^{M}\\ \vdots & \vdots & \cdots &\vdots \\ 1 & \mathrm x_N & \cdots &\mathrm x_N^M \end{matrix} \right\}\tag{5}$$

利用正规方程法，即公式(5)，我们可以得到如下$M$取不同值时的函数表达式：

图1

图2

图3

图4

---

图1-图4表明，随着$M$的增大，函数图像对训练样本的拟合越来越好，即训练误差越来越小。但是很明显图3的图像与原始的正弦函数图像最相似，即预测误差最小。 下表给出了图1-图4对应的最优的$\mathrm w$ 的取值：

w	$M=0$	$M=1$	$M=3$	$M=9$
$\omega_0$	$-0.0379$	$0.8309$	$-0.2655$	$-6.5887*10^{-2}$
$\omega_1$		$-1.9631$	$13.1817$	$-1.9234*10^1$
$\omega_2$			$-38.3154$	$5.2109*10^2$
$\omega_3$			$25.9214$	$-3.8321*10^3$
$\omega_4$				$1.3080*10^4$
$\omega_5$				$-2.1917*10^4$
$\omega_6$				$1.2754*10^4$
$\omega_7$				$1.1027*10^4$
$\omega_8$				$-1.8864*10^4$
$\omega_9$				$7.2725*10^3$

机器学习的目的就是选取最优的$M$值，最小化预测误差。但是实际值中，预测误差是在算法之后才能得到的(不然的话，预测有什么用)，我们都是通过验证误差来模拟预测误差。也就是说，我们一般把已经标记的数据集分为训练集和验证集，通过训练集来得到给定不同$M$时最小验证误差，从而选择最佳的$M$。图5和图6给出了$M$取值不同情况下的训练误差和验证误差：

图5 10个训练样例

图6 100个训练样例

在图5和图6中，训练样例和验证样例都是由公式(3)给出，但是图5只有10个训练样例，图6有100个训练样例，验证样例都为100。从图中可知，训练误差都是随$M$增加而下降。图5中，当训练样例为10个时，此时我们可以选择$M=3$或者$M=6$得到较小的验证误差。当训练样例足够多时，如图6所示，此时$M$越大，验证误差越好。

根据上述讨论，我们可以总结如下：

当训练样例较少时，我们需要选择合适的模型的复杂度，即这里的$M$值；当训练样例较多时，我们选择的模型越复杂越好，即选择较大的$M$值。

---

防止过拟合

当训练数据较少，而模型较为复杂时，容易出现过拟合。如在图1-图6中，只有$10$ 个训练数据，当$M=9$时，误差变大，这时出现过拟合现象。因此，我们可以通过增加训练数据和正则化来防止过拟合。

增加训练数据

图7和图8给出了，在$M=9$情况下，不同训练样例对函数表达式的模拟情况。可见，当训练样例较多时，得到的模型与原始模型(正选函数)更接近。

图7 50个训练样例

图8 100个训练样例

正则化

从上面的表格中可以看出，当过拟合时($M=9$)，输入变量$\mathrm x$的系数$\mathrm{w}$的系数变得特别大.此时，当$\mathrm{x}$变动十分小时，输出$\hat y$也变得很大，这就导致了预测时误差变大。此时，我们可以对误差函数加入惩罚项，来限制$\mathrm{w}$的取值：

$$E=\sum\limits_{i=1}^N{(\hat y_i-y_i)^2}+\frac{\lambda}{2}\lvert\mathrm{w}\rvert^2\tag{6}$$

公式(6)中的$\lambda$可以自己调节来选取合适的值。

同样地，我们可以使用正规方程来使得新的误差函数(6)最小。此时的解析解可以表示为：

$$\mathrm{\bar w}=[\bar\phi^{\mathrm T}(\mathrm X)\bar\phi(\mathrm X)+\lambda\mathrm I_0]^{-1}\bar\phi^{\mathrm T}(\mathrm X)\mathrm y\tag{7}$$

其中$\mathrm I_0$是一个$(M+1)\times(M+1)$的对角矩阵，且第一个对角元素为0(因为我们一般在正则项中不考虑$\omega_0$，见P10, PRML)，其他对角元素为$1$。

注意：公式(7)可以对$E$求导，使其为零得到，这里就不详细推导。

图9

图10

---

图9和图10给出了当测试样例为10，$\lambda$取不同值的函数拟合情况。可以看出，$\lambda$越小，对训练样例的拟合越好，即训练误差越小，但是此时图像与原始的正弦函数差别较大。当$\lambda=0$时，即不考虑正则化，此时对应的是图4。图9和图10说明，我们可以通过加入正则项来避免过拟合的情况。

---

偏差与方差

理论推导

机器学习的目的就是最小化误差。一般采用的误差，如线性回归的平方和误差，逻辑回归的交叉熵误差。这些误差都是假定训练样例的权重一样，但是实际中，每个样例出现的概率是不同的。因此，我们这里定义一个平均误差函数:

$$\mathbb E[E]=\int\int E(\hat y,y)p(\mathrm x,y)d\mathrm xdy\tag{8}$$

这里的$E(\hat y,y)$就是我们常用的误差函数,如下：

$$E(\hat y, y)=(\hat y_i-y_i)^2\tag{9}$$ $$E(\hat y, y)=-[y_i\log{\hat y_i}+(1-y_i)\log{(1-\hat y_i)}]\tag{10}$$

可见，公式(8)是一个广义的误差函数。这里我们平方和误差函数为例，将公式(9)带入公式(8)中，我们有

$$\mathbb E[E]=\int\int (\hat y-y)^2p(\mathrm x,y)d\mathrm xdy\tag{11}$$

通过求导，令导数为0，我们可以得到最佳的函数表达式：

$$\frac{\partial\mathbb E(E)}{\partial\hat y}=2\int(\hat y-y)p(\mathrm x,y)dy=0\tag{12}$$ $$\hat y=\frac{\int yp(\mathrm x,y)dy}{p(\mathrm x)}=\int yp(y\mid\mathrm x)dy=\mathbb E_y[y\mid\mathrm x]=y^\star\tag{13}$$

根据公式(13)我们可以重写(11)

$$\begin{align} \mathbb E(E)&=\int\int(\hat y-y^\star+y^\star-y)^2p(\mathrm x,y)d\mathrm xdy\nonumber\\ &=\int\int[(\hat y-y^\star)^2+2(\hat y-y^\star)(y^\star-y)+(y^\star-y)^2]d\mathrm xdy\nonumber\\ &=\int(\hat y-y^\star)^2p(\mathrm x)d\mathrm x+\int(y^\star-y)^2p(\mathrm x)d\mathrm x\tag{14} \end{align}$$

公式(14)中的第二项与我们要求的函数表达式$\hat y$没有关系。因此，当我们得到最优的函数表达式($\hat y=y^\star$)，即公式第一项为0，第二项即为我们得到的最小误差值。然而，由于训练数据有限(一般假定训练数据集为$\mathcal D={\mathrm x_i\mid i=1,2,\cdots,N}$)，得到最优解($\hat y=y^\star=\mathbb E(y\mid\mathrm x)$)一般是比较困难的。但我们有充足的训练数据$\mathrm x$, 我们理论上可以得到条件期望$\mathbb E(y\mid\mathrm x)=\int yp(y\mid\mathrm x)dy$, 也就是最优的函数表达式$\hat y$。

真实的$y$ 与$\mathrm x$的关系由$p(\mathrm x,y)$决定，假定由$p(\mathrm x,y)$产生很多不同训练数据集 $\mathcal D$ 。对于每一个数据集$\mathcal D$, 我们都能通过机器学习算法得到一个函数表达式$\hat y{\mathcal D}$。那么，我们需要在所有可能的训练数据集来评价$\hat y{\mathcal D}$的好坏，即我们需要计算$\hat y_{\mathcal D}$在所以训练集上的平均误差。那么公式(14)可以写成：

$$\begin{align} \mathbb E(E)&=\int\mathbb E_{\mathcal D}[(\hat y-y^\star)^2]p(\mathrm x)d\mathrm x+\int(y^\star-y)^2p(\mathrm x)d\mathrm x\nonumber\\ &=\int\{(\mathbb E_{\mathcal D}[\hat y_{\mathcal D}]-y^\star)^2+\mathbb E_{\mathcal D}[(\hat y_{\mathcal D}-\mathbb E_{\mathcal D}[\hat y_{\mathcal D}])^2]\}p(\mathrm x)d\mathrm x+\int(y^\star-y)^2p(\mathrm x)d\mathrm x\nonumber\\ &=\underbrace{\int(\mathbb E_{\mathcal D}[\hat y_{\mathcal D}]-y^\star)^2p(\mathrm x)d\mathrm x}_{(bias)^2}+\underbrace{\int\mathbb E_{\mathcal D}[(\hat y_{\mathcal D}-\mathbb E_{\mathcal D}[\hat y_{\mathcal D}])^2]p(\mathrm x)d\mathrm x}_{variance}+\underbrace{\int(y^\star-y)^2p(\mathrm x)d\mathrm x}_{noise}\tag{15} \end{align}$$

可见，误差由bias (偏差)，variance (方差)，和noise (噪声)三部分组成。其中，bias和variance都和我们的模型选择$\hat y{\mathcal D}$有关。第三项noise可以表示为$\int(\mathbb E{\mathcal D}[y\mid\mathrm x]-y)^2p(\mathrm x)d\mathrm x$, 即可以看成是训练数据自身的特征：$y$的方差。对于公式(15)，我们可以计算各部分值如下：

• $\mathbb E{\mathcal D}[\hat y{\mathcal D}]$指的是对$\hat y$在$L$个数据集上求平均值，那么我们有：

  $$\mathbb E_{\mathcal D}[\hat y_{\mathcal D}(\mathrm x)]=\frac{1}{L}\sum\limits_{l=1}^{L}\hat y_l(\mathrm x)\tag{16}$$

• 将公式(16)带入公式(15)中，我们有：

$$(bias)^2=\frac{1}{N}\sum\limits_{i=1}^{N}\{\mathbb E_{\mathcal D}[\hat y_{\mathcal D}(\mathrm x_i)]-y^\star(\mathrm x_i)\}^2\tag{17}$$ $$variance=\frac{1}{N}\sum\limits_{i=1}^{N}\frac{1}{L}\sum\limits_{l=1}^{L}\{\hat y_l(\mathrm x_i)-\mathbb E_{\mathcal D}[\hat y_{\mathcal D}(\mathrm x_i)]\}^2\tag{18}$$

注意：在计算公式(15)的积分项$\int p(\mathrm x)d\mathrm x$ 时，我们采用的是将所有的 $\mathrm x$ 所得到的结果加和求平均，即$\sum/N$。因为这里我们假定所有训练样例都是均匀采样的。

一般来说，模型($\hat y(\mathrm x)$)越复杂，偏差越小，方差越大。因为模型越复杂，对于我们训练样例集的每一个样例$\mathrm x_i$的拟合较好，也就是说$\hat y_l(\mathrm x_i)$与$y^\star(\mathrm x_i)$比较接近，即公式(17)的值较小(偏差较小)；而此时，不同训练样例集产生的$\hat y_l(\mathrm x_i)$之间值存在较大波动，即公式(18)的值较大(方差较大)。换句话说，公式(17)，即偏差，针对的是$\hat y_l(\mathrm x_i)$与真实函数$y^\star(\mathrm x_i)$之间的误差；而公式(18)，即方差，针对的是不同数据集所得到的函数$\hat y_l(\mathrm x_i)$之间的误差。

---

偏差与方差的折中关系

我们通过函数表达式(3)来产生$L=100$个训练数据集，每个数据集包含$N=30$个训练样例。那么此时，我们有$y^\star=\sin(2\pi\mathrm x)$且此时公式(15)的第三项noise可以由$\epsilon$的分布求出。我们假定用$M=9$的多项式函数(2)来作为$\hat y$的表达式。当正则项的系数分别为$\lambda=10^{-3},10^{-12}$时，我们可以分别得到图11-12，和图13-14：

图11

图12

图13

图14

图11和图13表示的是在这100个训练数据集下$\hat y$关于$\mathrm x$ 的函数图像。图12和图14表示的是在左图100条函数取平均情况下的函数图像，其中红色曲线是我们最优的函数$y^\star$。左图可以反映各个函数表达式间的差别，公式(18)， 即方差，右图表示的是预测函数$\mathbb E{\mathcal D}[\hat y{\mathcal D}]$与最优函数$y^\star$的差别，公式(17)，即偏差。通过图11-14，我们可以看出偏差与方差的折中关系。

---

经过上面的分析，我们可以看出误差主要由偏差、方差和噪声组成，并从中可以看出模型的选择(e.g.,这里$M,\lambda$的选择)对误差的本质影响，从而指导模型的选择。由公式(16)-(18)可以看出，误差的分析是建立在很多数据集上的统计平均值。但是在实际中，训练数据集很少。当我们有很多的数据集，我们可以把它们看出一个大的数据集，这样我们就可以防止过拟合现象(见图7-8)。

---

附录

下面我们给出图1-图14的python源代码。注意，在运行代码时，可以自行调整自变量$M,N,\lambda$等。

图1-4的python源代码：

# -*- coding: utf-8 -*-
# @Time : 2020/4/16 23:40
# @Author : tengweitw

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn import datasets


# Set the format of labels
def LabelFormat(plt):
    ax = plt.gca()
    plt.tick_params(labelsize=14)
    labels = ax.get_xticklabels() + ax.get_yticklabels()
    [label.set_fontname('Times New Roman') for label in labels]
    font = {'family': 'Times New Roman',
            'weight': 'normal',
            'size': 16,
            }
    return font


def Polynomial_regression_normal_equation(train_data, train_target, test_data, test_target):
    # the 1st column is 1 i.e., x_0=1
    X = np.ones([np.size(train_data, 0), 1])
    X_test = np.ones([np.size(test_data, 0), 1])
    # Here to change M !!!!!!!
    M = 2

    for i in range(1, M + 1):
        temp = train_data ** i
        temp_test = test_data ** i
        X = np.concatenate((X, temp), axis=1)
        X_test = np.concatenate((X_test, temp_test), axis=1)
    # X is a 10*M-dim matrix

    # Normal equation
    w_bar = np.matmul(np.linalg.pinv(np.matmul(X.T, X)), np.matmul(X.T, train_target))

    # Training Error
    y_predict_train = np.matmul(X, w_bar)
    E_train = np.linalg.norm(y_predict_train - train_target) / len(y_predict_train)

    # Predicting
    y_predict_test = np.matmul(X_test, w_bar)

    # Prediction Error
    E_test = np.linalg.norm(y_predict_test - test_target) / len(y_predict_test)

    return y_predict_test, E_train, E_test


if __name__ == '__main__':
    # keep the same random training data
    seed_num = 100
    np.random.seed(seed_num)
    # 10 training data
    train_data = np.random.uniform(0, 1, (10, 1))
    train_data = np.sort(train_data, axis=0)

    np.random.seed(seed_num)
    train_target = np.sin(2 * np.pi * train_data) + 0.1 * np.random.randn(10, 1)

    test_data = np.linspace(0, 1, 100).reshape(100, 1)
    np.random.seed(seed_num)
    test_target = np.sin(2 * np.pi * test_data) + 0.01 * np.random.randn(100, 1)

    y_predict_test, E_train, E_test = Polynomial_regression_normal_equation(train_data, train_target, test_data,
                                                                            test_target)

    plt.figure()
    plt.plot(train_data, train_target, 'ro')
    plt.plot(test_data, y_predict_test, 'b-')

    # Set the labels
    font = LabelFormat(plt)
    plt.xlabel('x', font)
    plt.ylabel('y', font)
    plt.legend(['Training target', 'Predicted target,M=2'])
    plt.ylim([-1, 1])
    plt.show()

图5-8的python源代码：

# -*- coding: utf-8 -*-
# @Time : 2020/4/18 11:56
# @Author : tengweitw

import numpy as np
import matplotlib.pyplot as plt



# Set the format of labels
def LabelFormat(plt):
    ax = plt.gca()
    plt.tick_params(labelsize=14)
    labels = ax.get_xticklabels() + ax.get_yticklabels()
    [label.set_fontname('Times New Roman') for label in labels]
    font = {'family': 'Times New Roman',
            'weight': 'normal',
            'size': 16,
            }
    return font


def Polynomial_regression_normal_equation(train_data, train_target, cv_data, cv_target,test_data,M):
    # the 1st column is 1 i.e., x_0=1
    X = np.ones([np.size(train_data, 0), 1])
    X_cv = np.ones([np.size(cv_data, 0), 1])
    X_test = np.ones([np.size(test_data, 0), 1])

    for i in range(1, M + 1):
        temp = train_data ** i
        temp_cv = cv_data ** i
        temp_test = test_data ** i
        X = np.concatenate((X, temp), axis=1)
        X_cv = np.concatenate((X_cv, temp_cv), axis=1)
        X_test = np.concatenate((X_test, temp_test), axis=1)
    # X is a 10*M-dim matrix

    # Normal equation
    w_bar = np.matmul(np.linalg.pinv(np.matmul(X.T, X)), np.matmul(X.T, train_target))

    # Training Error
    y_predict_train = np.matmul(X, w_bar)
    E_train = np.linalg.norm(y_predict_train - train_target) / len(y_predict_train)

    # cross validation
    y_predict_cv = np.matmul(X_cv, w_bar)

    # prediction
    y_predict_test=np.matmul(X_test, w_bar)

    # Prediction Error
    E_cv = np.linalg.norm(y_predict_cv - cv_target) / len(y_predict_cv)
    print(w_bar)
    return y_predict_test, y_predict_cv, E_train, E_cv


if __name__ == '__main__':
    # keep the same random training data
    seed_num = 100
    np.random.seed(seed_num)

    # training data
    num_training=50
    train_data = np.random.uniform(0, 1, (num_training, 1))
    train_data = np.sort(train_data, axis=0)
    np.random.seed(seed_num)
    train_target = np.sin(2 * np.pi * train_data) + 0.1 * np.random.randn(num_training, 1)

    # 100 cross validation data
    num_cv=100
    cv_data = np.random.uniform(0, 1, (num_cv, 1))
    cv_data = np.sort(cv_data, axis=0)
    np.random.seed(seed_num)
    cv_target = np.sin(2 * np.pi * cv_data) + 0.1 * np.random.randn(num_cv, 1)

    # testing data
    test_data = np.linspace(0, 1, 100).reshape(100, 1)

    M=9+1
    E_train=np.zeros((M,1))
    E_cv=np.zeros((M,1))
    # change M
    for i in range(M):
        y_predict_test,y_predict_cv, E_train[i], E_cv[i] = Polynomial_regression_normal_equation(train_data, train_target, cv_data, cv_target,test_data, i)

    plt.figure()
    plt.plot(E_train, 'r-o')
    plt.plot(E_cv,'b-s')


    # Set the labels
    font = LabelFormat(plt)
    plt.xlabel('$M$', font)
    plt.ylabel('Error', font)
    plt.legend(['Training Error', 'Cross Error'],loc='upper center')
    plt.grid()
    plt.show()

    plt.figure()
    plt.plot(train_data, train_target, 'ro')
    plt.plot(test_data, y_predict_test, 'b-')

    # Set the labels
    font = LabelFormat(plt)
    plt.xlabel('x', font)
    plt.ylabel('y', font)
    plt.legend(['Training target', 'Predicted target,M=9'])
    plt.ylim([-1, 1])
    plt.show()

图9-10的python源代码：

# -*- coding: utf-8 -*-
# @Time : 2020/4/18 23:12
# @Author : tengweitw


import numpy as np
import matplotlib.pyplot as plt



# Set the format of labels
def LabelFormat(plt):
    ax = plt.gca()
    plt.tick_params(labelsize=14)
    labels = ax.get_xticklabels() + ax.get_yticklabels()
    [label.set_fontname('Times New Roman') for label in labels]
    font = {'family': 'Times New Roman',
            'weight': 'normal',
            'size': 16,
            }
    return font


def Polynomial_regression_normal_equation(train_data, train_target, cv_data, cv_target,test_data,M):
    # the 1st column is 1 i.e., x_0=1
    X = np.ones([np.size(train_data, 0), 1])
    X_cv = np.ones([np.size(cv_data, 0), 1])
    X_test = np.ones([np.size(test_data, 0), 1])
    # Here to change lambda
    Lambda=1e-12
    I0= np.eye(M+1)
    I0[0]=0

    for i in range(1, M + 1):
        temp = train_data ** i
        temp_cv = cv_data ** i
        temp_test = test_data ** i
        X = np.concatenate((X, temp), axis=1)
        X_cv = np.concatenate((X_cv, temp_cv), axis=1)
        X_test = np.concatenate((X_test, temp_test), axis=1)
    # X is a 10*M-dim matrix

    # Normal equation
    w_bar = np.matmul(np.linalg.pinv(np.matmul(X.T, X)+Lambda*I0), np.matmul(X.T, train_target))

    # Training Error
    y_predict_train = np.matmul(X, w_bar)
    E_train = np.linalg.norm(y_predict_train - train_target) / len(y_predict_train)

    # cross validation
    y_predict_cv = np.matmul(X_cv, w_bar)

    # prediction
    y_predict_test=np.matmul(X_test, w_bar)

    # Prediction Error
    E_cv = np.linalg.norm(y_predict_cv - cv_target) / len(y_predict_cv)
    print(w_bar)
    return y_predict_test, y_predict_cv, E_train, E_cv


if __name__ == '__main__':
    # keep the same random training data
    seed_num = 100
    np.random.seed(seed_num)

    # training data
    num_training=10
    train_data = np.random.uniform(0, 1, (num_training, 1))
    train_data = np.sort(train_data, axis=0)
    np.random.seed(seed_num)
    train_target = np.sin(2 * np.pi * train_data) + 0.1 * np.random.randn(num_training, 1)

    # 100 cross validation data
    num_cv=100
    cv_data = np.random.uniform(0, 1, (num_cv, 1))
    cv_data = np.sort(cv_data, axis=0)
    np.random.seed(seed_num)
    cv_target = np.sin(2 * np.pi * cv_data) + 0.1 * np.random.randn(num_cv, 1)

    # testing data
    test_data = np.linspace(0, 1, 100).reshape(100, 1)

    M=9

    y_predict_test,y_predict_cv, E_train, E_cv = Polynomial_regression_normal_equation(train_data, train_target, cv_data, cv_target,test_data, M)


    plt.figure()
    plt.plot(train_data, train_target, 'ro')
    plt.plot(test_data, y_predict_test, 'b-')

    # Set the labels
    font = LabelFormat(plt)
    plt.xlabel('x', font)
    plt.ylabel('y', font)
    plt.legend(['Training target', 'Predicted target,$M=9,\lambda=10^{-12}$'])
    plt.ylim([-1, 1])
    plt.show()

图11-14的python源代码：

# -*- coding: utf-8 -*-
# @Time : 2020/4/20 10:46
# @Author : tengweitw


import numpy as np
import matplotlib.pyplot as plt



# Set the format of labels
def LabelFormat(plt):
    ax = plt.gca()
    plt.tick_params(labelsize=14)
    labels = ax.get_xticklabels() + ax.get_yticklabels()
    [label.set_fontname('Times New Roman') for label in labels]
    font = {'family': 'Times New Roman',
            'weight': 'normal',
            'size': 16,
            }
    return font


def Polynomial_regression_normal_equation(train_data, train_target,test_data,M):
    # the 1st column is 1 i.e., x_0=1
    X = np.ones([np.size(train_data, 0), 1])
    X_test = np.ones([np.size(test_data, 0), 1])
    # Here to change lambda
    Lambda=1e-8
    I0= np.eye(M+1)
    I0[0]=0

    for i in range(1, M + 1):
        temp = train_data ** i
        temp_test = test_data ** i
        X = np.concatenate((X, temp), axis=1)
        X_test = np.concatenate((X_test, temp_test), axis=1)
    # X is a 10*M-dim matrix

    # Normal equation
    w_bar = np.matmul(np.linalg.pinv(np.matmul(X.T, X)+Lambda*I0), np.matmul(X.T, train_target))

    # Training Error
    y_predict_train = np.matmul(X, w_bar)
    E_train = np.linalg.norm(y_predict_train - train_target) / len(y_predict_train)

    # prediction
    y_predict_test=np.matmul(X_test, w_bar)

    # Prediction Error

    return y_predict_test, E_train


if __name__ == '__main__':
    # L is  the number of training data sets
    L=100

    # number of each training data
    N=30


    plt.figure()
    # For L training datasets
    for l in range(L):

        train_data = np.random.uniform(0, 1, (N, 1))
        train_data = np.sort(train_data, axis=0)
        train_target = np.sin(2 * np.pi * train_data) + 0.1* np.random.randn(N, 1)
        # testing data
        test_data = np.linspace(0, 1, 100).reshape(100, 1)

        M=9
        y_predict_test, E_train = Polynomial_regression_normal_equation(train_data, train_target, test_data, M)

        if l==0:
            predict_target=y_predict_test
        else:
            predict_target=np.hstack((y_predict_test,predict_target))

        # plt.plot(train_data, train_target, 'ro')
        plt.plot(test_data, y_predict_test, 'b-')

    # Set the labels
    font = LabelFormat(plt)
    plt.xlabel('x', font)
    plt.ylabel('y', font)
    plt.legend([ 'Predicted target,$M=9,\lambda=10^{-3}$'])

    plt.ylim(-1,1)
    plt.show()

    predict_target_avg=np.mean(predict_target,axis=1)
    plt.figure()
    plt.plot(test_data,predict_target_avg,'r-')
    test_target = np.sin(2 * np.pi * test_data)
    plt.plot(test_data,test_target,'b-')

    # Set the labels
    font = LabelFormat(plt)
    plt.xlabel('x', font)
    plt.ylabel('y', font)
    plt.legend(['True model', 'Predicted model,$M=9,\lambda=10^{-3}$'])

    plt.show()