机器学习——多项式拟合

文章目录

一、多项式拟合
- 1.1 多项式函数拟合
- 1.2 过拟合与正则化
- 1.3 维数灾难
二、基于概率的曲线拟合
- 2.1 高斯分布
- 2.2 重新理解曲线拟合
- - 最大似然
  - 最大化后验概率（maximum posterior, MAP）
  - Bayes 曲线拟合

开发环境jupyter notebook

一、多项式拟合

引入：通过这个示例，从最简单的机器学习的算法入手实现复现流程

假设我们有两个实值变量 x,tx,tx,t，满足关系：t=sin(2πx)+ϵt=sin(2πx)+ϵt=sin(2πx)+ϵ，其中 ϵ 是一个服从高斯分布的随机值。
假设我们有 NNN 组 (x,t)(x,t)(x,t) 的观测值 x≡(x1,…,xN)T,t≡(t1,…,tN)Tx≡(x_1,…,x_N)^T,t≡(t_1,…,t_N)^Tx≡(x1,…,xN)T,t≡(t1,…,tN)T：

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
%matplotlib inlineN = 10                                                        # 设置 N
x_tr = np.linspace(0, 1, N)                                   # 生成 0，1 之间等距的 N 个 数
t_tr = np.sin(2 * np.pi * x_tr) + 0.25 * np.random.randn(N)   # 计算 t,引入随机噪声# 绘图
xx = np.linspace(0, 1, 500)fig, ax = plt.subplots()
ax.plot(x_tr, t_tr, 'co')
ax.plot(xx, np.sin(2 * np.pi * xx), 'g')
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-1.5, 1.5)
ax.set_xticks([0, 1])
ax.set_yticks([-1, 0, 1])
ax.set_xlabel("$x$", fontsize="x-large")
ax.set_ylabel("$t$", fontsize="x-large")plt.show()

使用这 NNN 个数据点作为训练集，我们希望得到这个一个模型：给定一个新的输入x^\hat{x}x^，预测他对应的输出 t^\hat{t}t^。

1.1 多项式函数拟合

我们使用曲线拟合的方法来拟合这样一个多项式函数：
y(x,w)=w0+w1x+w0+w2x2+⋯+wMxM=∑j=0Mwjxj(1)y(\mathbf{x},\mathbf{w})=w_0+w_1x+w_0+w_2x^2 +\cdots +w_Mx^M = \sum_{j=0}^{M}w_jx^j \tag1y(x,w)=w0+w1x+w0+w2x2+⋯+wMxM=j=0∑Mwjxj(1)

其中 MMM 是多项式的阶数，xjx^jxj 表示 xxx 的 jjj 次方，w≡(w0,w1,…,wM)w≡(w_0,w_1,…,w_M)w≡(w0,w1,…,wM) 表示多项式的系数。

这些多项式的系数可以通过我们的数据拟合得到，即在训练集上最小化一个关于 y(x,w)y(x,w)y(x,w) 和 ttt 的损失函数。常见的一个损失函数是平方误差和，定义为：
E(w)=12∑i=1N{y(x,w)−tn}2(2)E(\mathbf{w}) = \frac{1}{2}\sum_{i=1}^{N}\left \{ y(x,\mathbf{w})-t_n \right \}^2 \tag2E(w)=21i=1∑N{y(x,w)−tn}2(2)

因子 12\frac{1}{2}21 是为了之后的计算方便加上的。这个损失函数是非负的，当且仅当函数 y(x,w)y(x,\mathbf{w})y(x,w) 通过所有的数据点时才会为零。
对于这个损失函数，因为它是一个关于 w\mathbf{w}w 的二次函数，其关于 w\mathbf{w}w 的梯度是一个关于 www 的线性函数，因此我们可以找到一个唯一解 w⋆\mathbf{w}^⋆w⋆。
∂E(w)wj=∑n=1N(∑j=0Mwjxnj−tn)xnj=0(3)\frac{\partial E(\mathbf{w})}{w_{j}} = \sum_{n=1}^{N}\left ( \sum_{j=0}^{M}w_jx_n^j - t_n \right )x_n^j = 0 \tag3wj∂E(w)=n=1∑N(j=0∑Mwjxnj−tn)xnj=0(3)

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
%matplotlib inlineN = 10                                                        # 设置 N
x_tr = np.linspace(0, 1, N)                                   # 生成 0，1 之间等距的 N 个 数
t_tr = np.sin(2 * np.pi * x_tr) + 0.25 * np.random.randn(N)   # 计算 t,引入随机噪声fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes = axes.flatten()
Ms = [0, 1, 3, 9]                                             # 拟合参数是多项式的阶数 M，我们可以看看当 M=0,1,3,9 时的效果for ax, M in zip(axes, Ms):coeff = np.polyfit(x_tr, t_tr, M)                         # 计算参数f = np.poly1d(coeff)                                      # 生成函数 y(x, w)# 绘图xx = np.linspace(0, 1, 500)ax.plot(x_tr, t_tr, 'co')ax.plot(xx, np.sin(2 * np.pi * xx), 'g')ax.plot(xx, f(xx), 'r')ax.set_xlim(-0.1, 1.1)ax.set_ylim(-1.5, 1.5)ax.set_xticks([0, 1])ax.set_yticks([-1, 0, 1])ax.set_xlabel("$x$",fontsize="x-large")ax.set_ylabel("$t$",fontsize="x-large")ax.text(0.6, 1, '$M={}$'.format(M), fontsize="x-large")
plt.show()

可以看到，M=3 似乎是其中较好的一个选择，M=9 虽然拟合的效果最好（通过了所有的训练数据点），但很明显过拟合了。

检测模型的效果，用一组与训练集相同分布的数据进行测试，然后计算不同模型选择下，训练集和测试集上的 E(w⋆)E(\mathbf{w}^⋆)E(w⋆) 值。

注意到随着测试点数目的变化，E(w⋆)E(\mathbf{w}^⋆)E(w⋆) 的尺度也在不断变化，因此一个更好的选择是使用 root-mean-square (RMS) 误差：ERMS=2E(wsstar)/NE_{RMS}=\sqrt{2E(\mathbf{w}^sstar)/N}ERMS=2E(wsstar)/N
RMS 误差的衡量尺度和单位与目标值 t 一致。

我们用相同的方法产生100组数据作为测试集，计算不同 M 下的 RMS 误差：

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
%matplotlib inlineN = 10                                                        # 设置 N
x_tr = np.linspace(0, 1, N)                                   # 生成 0，1 之间等距的 N 个 数
t_tr = np.sin(2 * np.pi * x_tr) + 0.25 * np.random.randn(N)   # 计算 t,引入随机噪声x_te = np.random.rand(100)
t_te = np.sin(2 * np.pi * x_te) + 0.25 * np.random.randn(100)
rms_tr, rms_te = [], []for M in range(10):coeff = np.polyfit(x_tr, t_tr, M)         # 计算参数f = np.poly1d(coeff)                      # 生成函数 y(x, w)# RMSrms_tr.append(np.sqrt(((f(x_tr) - t_tr) ** 2).sum() / x_tr.shape[0]))rms_te.append(np.sqrt(((f(x_te) - t_te) ** 2).sum() / x_te.shape[0]))# 画图
fig, ax = plt.subplots()ax.plot(range(10), rms_tr, 'bo-', range(10), rms_te, 'ro-')
ax.set_xlim(-1, 10)
ax.set_ylim(0, 1)
ax.set_xticks(range(0, 10, 3))
ax.set_yticks([0, 0.5, 1])
ax.set_xlabel("$M$",fontsize="x-large")
ax.set_ylabel("$E_{RMS}$",fontsize="x-large")
ax.legend(['Traning', 'Test'], loc="best")plt.show()

可以看到 M=9 时，虽然训练集上的误差已经降到 0，但是测试集上的误差却很大，为了更好的拟合这些点，系数都会变得很大：

for i, w in enumerate(np.polyfit(x_tr, t_tr, 9)):print "w_{}, {:.2f}".format(9 - i, w)'''
w_9, -69435.47
w_8, 322877.34
w_7, -631263.51
w_6, 673357.39
w_5, -424976.96
w_4, 160681.72
w_3, -34984.60
w_2, 3903.04
w_1, -158.55
w_0, -0.15
'''

1.2 过拟合与正则化

过拟合问题：当训练数据量 N 变多时，模型拟合的过拟合现象在减少。
当模型的复杂度固定时，随着数据的增加，过拟合的现象也在逐渐减少。
如下：M=9 的模型的表现：

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
%matplotlib inlinefig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes = axes.flatten()
M = 9for ax, N in zip(axes, (15, 150)):x_tr_more = np.linspace(0, 1, N)                                       # 生成 0，1 之间等距的N个数t_tr_more = np.sin(2 * np.pi * x_tr_more) + 0.25 * np.random.randn(N)  # 计算 tcoeff = np.polyfit(x_tr_more, t_tr_more, M)                            # 计算参数f = np.poly1d(coeff)                                                   # 生成函数 y(x, w)# 绘图xx = np.linspace(0, 1, 500)ax.plot(x_tr_more, t_tr_more, 'co')ax.plot(xx, np.sin(2 * np.pi * xx), 'g')ax.plot(xx, f(xx), 'r')ax.set_xlim(-0.1, 1.1)ax.set_ylim(-1.5, 1.5)ax.set_xticks([0, 1])ax.set_yticks([-1, 0, 1])ax.set_xlabel("$x$", fontsize="x-large")ax.set_ylabel("$t$", fontsize="x-large")ax.text(0.6, 1, '$N={}$'.format(N), fontsize="x-large")plt.show()

正则化：如果我们一定要在 N=10 的数据上使用 M=9 的模型，那么一个通常的做法是给参数加一个正则项的约束防止过拟合，一个最常用的正则项是平方正则项，即控制所有参数的平方和大小：
E(w)^=12∑i=1N{y(xn,w)−tn}2+λ2∣∣w∣∣2(4)\hat{E(\mathbf{w})} = \frac{1}{2}\sum_{i=1}^{N}\left \{ y(x_n,\mathbf{w})-t_n \right \}^2+\frac{\lambda }{2}||\mathbf{w}||^2 \tag4E(w)^=21i=1∑N{y(xn,w)−tn}2+2λ∣∣w∣∣2(4)

其中 ∥w∥2≡ww⊤w=w02+…∣wM2，λ∥\mathbf{w}∥^2≡w\mathbf{w}^⊤\mathbf{w}=w_0^2+…|w_M^2，λ∥w∥2≡ww⊤w=w02+…∣wM2，λ 是控制正则项和误差项的相对重要性。
若设向量 ϕ(x)ϕ(x)ϕ(x) 满足 ϕi(x)=xi,i=0,1,…,Mϕ_i(x)=x^i,i=0,1,…,Mϕi(x)=xi,i=0,1,…,M，则对 w\mathbf{w}w 最小化，解应当满足：
[∑n=1Nϕ(xn)ϕ(xn)……T+λI]w=∑n=1Ntnϕ(xn)(5)\left [ \sum_{n=1}^{N}\phi (x_n)\phi (x_n)……T+\lambda I \right ]\mathbf{w} = \sum_{n=1}^{N}t_n\phi (x_n) \tag5[n=1∑Nϕ(xn)ϕ(xn)……T+λI]w=n=1∑Ntnϕ(xn)(5)

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
%matplotlib inlinedef phi(x, M):return x[:,None] ** np.arange(M + 1)N = 10                                                        # 设置 N
x_tr = np.linspace(0, 1, N)                                   # 生成 0，1 之间等距的 N 个 数
t_tr = np.sin(2 * np.pi * x_tr) + 0.25 * np.random.randn(N)   # 计算 t,引入随机噪声# 加正则项的解
M = 9
lam = 0.0001phi_x_tr = phi(x_tr, M)
S_0 = phi_x_tr.T.dot(phi_x_tr) + lam * np.eye(M+1)
y_0 = t_tr.dot(phi_x_tr)coeff = np.linalg.solve(S_0, y_0)[::-1]
f = np.poly1d(coeff)                                          # 生成函数 y(x, w)xx = np.linspace(0, 1, 500)fig, ax = plt.subplots()
ax.plot(x_tr, t_tr, 'co')
ax.plot(xx, np.sin(2 * np.pi * xx), 'g')
ax.plot(xx, f(xx), 'r')
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-1.5, 1.5)
ax.set_xticks([0, 1])
ax.set_yticks([-1, 0, 1])
ax.set_xlabel("$x$", fontsize="x-large")
ax.set_ylabel("$t$", fontsize="x-large")plt.show()

1.3 维数灾难

在上述多项式拟合问题中只考虑了 xxx 是一维变量的结果，但事实上 xxx 的维度可能远远不止一维，考虑一个 DDD 维的输入 x\bf xx 并且选择阶数 M=3\bf M=3M=3，那么我们有
y(x,w)=w0+∑i=1Dwixi+∑i=1D∑j=1Dwijxixj+∑i=1D∑j=1D∑k=1Dwijxixjxk(6)y(\mathbf{x},\mathbf{w})=w_0+\sum_{i=1}^{D}w_ix_i+\sum_{i=1}^{D}\sum_{j=1}^{D}w_{ij}x_ix_j+\sum_{i=1}^{D}\sum_{j=1}^{D}\sum_{k=1}^{D}w_{ij}x_ix_jx_k \tag6y(x,w)=w0+i=1∑Dwixi+i=1∑Dj=1∑Dwijxixj+i=1∑Dj=1∑Dk=1∑Dwijxixjxk(6)

可以看到，随着 DDD 的增大，独立参数的个数也增大到 D3D^3D3 的量级。对于 MMM 阶的多项式，系数的量级将变成 DMD^MDM。
也就是说，随着 DDD 的增大，独立参数的数目呈幂数增长。

二、基于概率的曲线拟合

2.1 高斯分布

高斯分布（Gaussian distribution），又叫正态分布（normal distribution）。

对于实值变量 xxx，高斯分布定义为：
N(x∣μ,σ2)=1(2πσ2)exp⁡{−12σ2(x−μ)2}(7)\mathcal{N}\left(x\left|~\mu,\sigma^2\right.\right) = \frac{1}{\sqrt{(2\pi\sigma^2)}} \exp\left\{-\frac{1}{2\sigma^2}(x-\mu)^2\right\} \tag7N(x∣∣ μ,σ2)=(2πσ2)1exp{−2σ21(x−μ)2}(7)

参数为均值 μμμ 和方差 σ2σ^2σ2。方差的平方根 σσσ 叫做标准差，方差的倒数 β=1σ2β=\frac{1}{σ^2}β=σ21 叫做精度。

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.stats import norm
%matplotlib inlinexx = np.linspace(-3, 3, 200)norm_xx = norm.pdf(xx)fig, ax = plt.subplots()ax.plot(xx, norm_xx, "r")
ax.set_ylim(0, 0.5)
ax.set_ylabel(r"$\mathcal{N}\left(x|\mu,\sigma^2\right)$", fontsize="xx-large")
ax.set_yticks([])
ax.set_yticklabels([])ax.set_xticks([0])
ax.set_xticklabels([r"$\mu$"], fontsize="xx-large")ax.text(-.1, 0.25, "$2\sigma$", fontsize="xx-large")ax.annotate("",xy=(-1, 0.24), xycoords='data',xytext=(1, 0.24), textcoords='data',arrowprops=dict(arrowstyle="<->",connectionstyle="arc3"), )plt.show()

多维高斯分布
对于 D 维的向量 x，高斯分布定义为：
N(x∣μ,Σ)=1(2π)D1∣Σ∣exp⁡{−12(x−μ)⊤Σ−1(x−μ)}(8)\mathcal{N}\left(\mathbf x\left|~\mathbf{\mu, \Sigma}\right.\right) = \frac{1}{\sqrt[D]{(2\pi)}} \frac{1}{\sqrt{|\mathbf\Sigma|}} \exp \left\{-\frac{1}{2}(\mathbf x - \mathbf \mu)^\top\mathbf\Sigma^{-1}(\mathbf x - \mathbf \mu)\right\} \tag8N(x∣ μ,Σ)=D(2π)1∣Σ∣1exp{−21(x−μ)⊤Σ−1(x−μ)}(8)

其中，DDD 维向量 μμμ 是均值，D×DD×DD×D 矩阵 Σ\mathbf{Σ}Σ 是方差，∣Σ∣|\mathbf{Σ}|∣Σ∣ 是其行列式。

最大似然估计

假设我们现在有 NNN 组对 xxx 的观测数据 x=(x1,…,xN)T\mathsf x = (x_1,\dots,x_N)^{\text T}x=(x1,…,xN)T，这些数据是独立同分布（independent and identically distributed, i.i.d.）的，都服从一个均值 μμμ，方差 σ2σ^2σ2 的高斯分布。那么在给定这些参数的情况下，出现这些观测数据的概率，或者从参数的角度来说，似然函数为：详细求解流程
p(x∣μ,σ2)=∏n=1NN(xn∣μ,σ2)(9)p(\mathsf x~|~\mu, \sigma^2)=\prod_{n=1}^N \mathcal N\left(x_n \left|~\mu, \sigma^2\right.\right) \tag9p(x ∣ μ,σ2)=n=1∏NN(xn∣∣ μ,σ2)(9)

符合高斯分布最大似然解：

对 μμμ 最大化（即样本均值）μ=1N∑n=1Nxn\mu = \frac 1 N \sum_{n=1}^N x_nμ=N1∑n=1Nxn
对 σ2σ^2σ2 最大化（即样本方差）：σ2=1N∑n=1N(xn−μ)2\sigma^2 = \frac 1 N \sum_{n=1}^N (x_n-\mu)^2σ2=N1∑n=1N(xn−μ)2

方差的一个无偏估计为：（显然：随着 N 的增大，方差估计的误差也随之增大。）
σ~2=NN−1σ2=1N−1∑n=1N(xn−μ)2(10)\tilde \sigma^2 = \frac{N}{N-1}\sigma^2=\frac{1}{N-1}\sum_{n=1}^N(x_n-\mu)^2 \tag{10}σ~2=N−1Nσ2=N−11n=1∑N(xn−μ)2(10)

2.2 重新理解曲线拟合

对于曲线拟合的问题，设训练集输入为x=(x1,…,xN)⊤\mathsf x=(x_1, \dots, x_N)^\topx=(x1,…,xN)⊤，对应的目标值为 t=(t1,…,tN)⊤\mathsf t=(t_1, \dots, t_N)^\topt=(t1,…,tN)⊤。

我们将我们的不确定性用高斯分布来表示，假设给定 xxx，对应的目标值 ttt 服从一个均值为 y(x,w)y(x,\mathbf w)y(x,w) 的高斯分布：
p(t∣x,w,β)=N(t∣y(x,w),β−1)(11)p(t\left|~x,\mathbf w,\beta\right.)=\mathcal N\left(t\left|~y(x,\mathbf w), \beta^{-1}\right.\right) \tag{11}p(t∣ x,w,β)=N(t∣∣ y(x,w),β−1)(11)

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.stats import norm
%matplotlib inlinexx = np.linspace(-0.9, 0.9, 100)
yy = 4 * xx - np.sin(xx * np.pi)fig, ax = plt.subplots()
ax.plot(xx, yy, color="red")ax.set_xlim(-1, 1)
ax.set_ylim(-4, 4)ax.set_xticks([0])
ax.set_xticklabels([r'$x_0$'], fontsize="xx-large")
ax.set_yticks([0])
ax.set_yticklabels([r'$y(x_0, \mathbf{w})$'], fontsize="xx-large")xx = np.linspace(-4, 4, 100)
yy = norm.pdf(xx, scale=0.5) / 5ax.plot([-1, 0], [0, 0], "g--")
ax.plot([0, 0], [-4, 4], "k")
ax.plot(yy, xx)ax.annotate("",xy=(0.75, -0.5), xycoords='data',xytext=(0.75, 0.5), textcoords='data',arrowprops=dict(arrowstyle="<->",connectionstyle="arc3"), )ax.text(0.77, -0.2, r'$2\sigma$', fontsize="xx-large")
ax.text(0.15, -1, r'$p(t|x_0,\mathbf{w}, \beta)$', fontsize="xx-large")
ax.text(0.5, 3, r'$y(x, \mathbf{w})$', fontsize="xx-large")plt.show()

最大似然

设训练集数据是独立同分布的，那么似然函数为：
p(t∣x,w,β)=∑i=1NN(tn∣y(x,w),β−1)(12)p(\mathsf t \left|\mathsf x, \mathbf w, \beta\right.) = \sum_{i=1}^N \mathcal{N}\left(t_n\left|~y(x,\mathbf w), \beta^{-1}\right.\right) \tag{12}p(t∣x,w,β)=i=1∑NN(tn∣∣ y(x,w),β−1)(12)

对数似然为：
ln⁡p(t∣x,w,β)=−β2∑n=1N{y(x,w)−tn}2+N2ln⁡β−N2ln⁡(2π)(13)\ln p(\mathsf t \left|\mathsf x, \mathbf w, \beta\right.) = -\frac{\beta}{2} \sum_{n=1}^N \{y(x,\mathbf w) - t_n\}^2 + \frac N 2 \ln \beta - \frac N 2 \ln(2\pi) \tag{13}lnp(t∣x,w,β)=−2βn=1∑N{y(x,w)−tn}2+2Nlnβ−2Nln(2π)(13)

设系数的最大似然解为 w\mathbf{w}w，从最大化对数似然的角度来看，求它的问题相当于最小化：
12∑n=1N{y(x,w)−tn}2(14)\frac{1}{2} \sum_{n=1}^N \{y(x,\mathbf w) - t_n\}^2 \tag{14}21n=1∑N{y(x,w)−tn}2(14)

这就是之前最小化平方误差和的结果。因此最小化平方误差和可以看成是高斯噪音假设下的最大似然的结果。
再对精度 βββ 求最大似然，我们有（可以理解为照搬之前求 σ2σ^2σ2 的结果）：
1βML=1N∑i=1N{y(x,w)−tn}2(15)\frac{1}{\beta_{ML}} = \frac{1}{N}\sum_{i=1}^N \{y(x,\mathbf w) - t_n\}^2 \tag{15}βML1=N1i=1∑N{y(x,w)−tn}2(15)

我们有了最大似然的结果之后，对于一个新的输入 xxx，其输出 ttt 应当满足：
p(t∣x,w,β)=N(t∣y(x,w),β−1)(16)p\left(t\left|~x,\mathbf w, \beta\right.\right) = \mathcal N\left(t\left|~y(x,\mathbf w), \beta^{-1}\right.\right) \tag{16}p(t∣ x,w,β)=N(t∣∣ y(x,w),β−1)(16)

最大化后验概率（maximum posterior, MAP）

假设我们对系数 w\mathbf ww 有一个先验的知识（MMM 是多项式阶数，加上常数项一共 M+1M+1M+1 维）：
p(w∣α)=N(w∣0,α−1I)=(α2π)(M+1)/2exp⁡{−α2w⊤w}(17)p(\mathbf w~|~\alpha) = \mathcal{N}(\mathbf w~|~\mathbf 0, \alpha^{-1} I) = \left(\frac{\alpha}{2\pi}\right)^{(M+1)/2} \exp\left\{-\frac{\alpha}{2}\mathbf{w}^\top\mathbf{w}\right\} \tag{17}p(w ∣ α)=N(w ∣ 0,α−1I)=(2πα)(M+1)/2exp{−2αw⊤w}(17)
ααα 控制这个模型的先验分布，这一类的参数通常被叫做超参（hyperparameters），Bayes 公式告诉我们，后验概率正比与先验概率和似然函数的乘积：
p(w∣x,t,α,β)∝p(t∣x,w,β)p(w∣α)(18)p(\mathbf w~|~\mathsf{x, t}, \alpha, \beta) \propto p(\mathsf t~|~\mathsf x, \mathbf w, \beta)~p(\mathbf w~|~\alpha) \tag{18}p(w ∣ x,t,α,β)∝p(t ∣ x,w,β) p(w ∣ α)(18)

我们可以通过最大化后验概率（maximum posterior, MAP）来决定参数 w\mathbf ww 的值，对上式求对数，并去掉跟 w\mathbf ww 无关的项，我们相当于要最大化：
−β2∑n=1N{y(x,w)−tn}2−α2w⊤w(19)-\frac{\beta}{2} \sum_{n=1}^N \{y(x,\mathbf w) - t_n\}^2 -\frac{\alpha}{2}\mathbf{w}^\top\mathbf{w} \tag{19}−2βn=1∑N{y(x,w)−tn}2−2αw⊤w(19)
即最小化
β2∑n=1N{y(x,w)−tn}2+α2w⊤w(20)\frac{\beta}{2} \sum_{n=1}^N \{y(x,\mathbf w) - t_n\}^2 +\frac{\alpha}{2}\mathbf{w}^\top\mathbf{w} \tag{20}2βn=1∑N{y(x,w)−tn}2+2αw⊤w(20)

因此，MAP 的结果相当于给多项式拟合加二范数正则化的结果，其中正则参数 λ=α/βλ=α/βλ=α/β。

Bayes 曲线拟合

虽然在 MAP 中，我们引入了先验分布，但是本质上它还是一个点估计，本质上并不是一个完全的 Bayes 估计。

一个完全的 Bayes 估计要求我们对 w\mathbf ww 的所有值进行积分。

对于之前的曲线拟合问题，给定训练集 xxx 和 ttt，对于一个新的测试样例 xxx，其目标值为 ttt，我们考虑预测的分布 p(t∣x,x,t)p(t | x,x,t)p(t∣x,x,t)（这里我们假定 β,αβ,αβ,α 两个参数已经给定了）

Bayes 公式给出：
p(t∣x,x,t)=∫p(t∣x,w)p(w∣x,t)dw(21)p(t~|~x,\mathsf{x, t}) = \int p(t~|~x,\mathbf w) p(w~|~\mathsf{x,t})d\mathbf w \tag{21}p(t ∣ x,x,t)=∫p(t ∣ x,w)p(w ∣ x,t)dw(21)

其中 p(t∣x,w)p(t~|~x,\mathbf w)p(t ∣ x,w) 是之前给定的高斯分布，p(w∣x,t)p(w~|~\mathsf{x,t})p(w ∣ x,t) 是训练集上的后验概率（也是一个高斯分布）。
由于高斯分布的性质，上面的式子本质上也是一个高斯分布，因此可以写成:
p(t∣x,x,t)=N(t∣m(x),s2(x))(22)p(t~|~x,\mathsf{x, t})=\mathcal{N}\left(t~|~m(x), s^2(x)\right) \tag{22}p(t ∣ x,x,t)=N(t ∣ m(x),s2(x))(22)

其中均值和方差分别为
m(x)=βϕ(x)TS∑n=1Nϕ(xn)tns2(x)=β−1+ϕ(x)⊤Sϕ(x)(23)\begin{aligned} m(x) & = \beta \phi(x)^\text{T} \mathbf S \sum_{n=1}^N \phi(x_n) t_n \\ s^2(x) & = \beta^{-1} + \phi(x)^\top \mathbf S \phi(x) \tag{23} \end{aligned} m(x)s2(x)=βϕ(x)TSn=1∑Nϕ(xn)tn=β−1+ϕ(x)⊤Sϕ(x)(23)

其中，ϕi(x)=xi,i=0,…,M\phi_i(x) = x^i, i = 0, \dots, Mϕi(x)=xi,i=0,…,M，矩阵 S：
S−1=αI+β∑n=1Nϕ(xn)⊤ϕ(xn)(24)\mathbf S^{-1}=\alpha I +\beta \sum_{n=1}^N \phi(x_n)^\top\phi(x_n) \tag{24}S−1=αI+βn=1∑Nϕ(xn)⊤ϕ(xn)(24)

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.stats import norm
%matplotlib inlinedef phi(x, M):return x[:,None] ** np.arange(M + 1)N = 10
x_tr = np.linspace(0, 1, N)                                   # 生成 0，1 之间等距的 N 个 数
t_tr = np.sin(2 * np.pi * x_tr) + 0.25 * np.random.randn(N)  # 计算 t# 加正则项的解
M = 9
alpha = 5e-3
beta = 11.1
lam = alpha / betaphi_x_tr = phi(x_tr, M)
A_0 = phi_x_tr.T.dot(phi_x_tr) + lam * np.eye(M+1)
y_0 = t_tr.dot(phi_x_tr)# 求解 Aw=y
coeff = np.linalg.solve(A_0, y_0)[::-1]
f = np.poly1d(coeff)# 绘图
xx = np.linspace(0, 1, 500)# Bayes估计的均值和标准差
S = np.linalg.inv(A_0 * beta)m_xx = beta * phi(xx, M).dot(S).dot(y_0)
s_xx = np.sqrt(1 / beta + phi(xx, M).dot(S).dot(phi(xx, M).T).diagonal())fig, ax = plt.subplots()
ax.plot(x_tr, t_tr, 'co')
ax.plot(xx, np.sin(2 * np.pi * xx), 'g')
ax.plot(xx, f(xx), 'r')
ax.fill_between(xx, m_xx-s_xx, m_xx+s_xx, color="pink")
ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-1.5, 1.5)
ax.set_xticks([0, 1])
ax.set_yticks([-1, 0, 1])
ax.set_xlabel("$x$", fontsize="x-large")
ax.set_ylabel("$t$", fontsize="x-large")plt.show()

下图粉红色部分就是 Bayes 估计给出的结果，红色曲线是 MAP 给出的结果。