线性模型参数求解的最大似然估计、MAP估计、正则最小二乘估计

\qquad 本文主要描述针对线性回归模型 $y(\mathbf{x},\mathbf{w})={\mathbf{w}}^T\mathbf{\varphi }(\mathbf{x})$ 采用最大似然估计、最大后验估计，以及正则最小二乘估计在求取参数 $\mathbf{w}$ 值的过程中的联系。
\qquad

1. 线性回归的概率模型

\qquad 广义的线性回归函数可定义为：

$y(\mathbf{x},\mathbf{w})=\sum _{j=1}^M{w}_j{\varphi }_j(\mathbf{x})={\mathbf{w}}^T\mathbf{\varphi }(\mathbf{x})$

\qquad 其中，权向量 $\mathbf{w}=[w_1,\cdots ,w_M{]}^T$，观测数据 $\mathbf{x}=[x_1,\cdots ,x_D{]}^T$，基函数 $\mathbf{\varphi }(\mathbf{x})=[{\varphi }_1(\mathbf{x}),\cdots ,{\varphi }_M(\mathbf{x}){]}^T$。

\qquad
\qquad 假设由线性回归函数 $y(\mathbf{x},\mathbf{w})$ 所刻画的目标变量 $t$ 满足：

$t=y(\mathbf{x},\mathbf{w})+\epsilon$

\qquad 已知观测误差 $\epsilon \sim \mathcal{N}(0,{\sigma }^2)$，即：$p(\epsilon )=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\epsilon }^2}{2{\sigma }^2}\right)$

\qquad

• 考虑一个观测数据集 $\mathbf{X}=[{\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N]$，对应目标值向量为 $\mathbf{t}=[t_1,\cdots ,t_N{]}^T$，那么第 $i$ 个目标变量 $t_i$ 满足：

$t_i=\sum _{j=1}^M{w}_j{\varphi }_j({\mathbf{x}}_i)+{\epsilon }_i$　或者　$t_i={\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)+{\epsilon }_i$

\qquad 其中，权向量 $\mathbf{w}=[w_1,\cdots ,w_M{]}^T$，第 $i$ 个观测数据 ${\mathbf{x}}_i=[x_{i1},\cdots ,x_{iD}{]}^T$。
\qquad

对于基本线性模型 $t=\sum _{j=1}^M{w}_j{x}_j+\epsilon$，此时 $M=D$
基函数满足 ${\varphi }_j(\mathbf{x})=x_j$，因而 $\mathbf{\varphi }({\mathbf{x}}_i)=[x_{i1},\cdots ,x_{ij},\cdots ,x_{iM}{]}^T={\mathbf{x}}_i$，那么：
　
$t_i=\sum _{j=1}^M{w}_j{x}_{ij}+{\epsilon }_i$　或者　$t_i={\mathbf{w}}^T{\mathbf{x}}_i+{\epsilon }_i$

\qquad

• 由于观测误差 ${\epsilon }_i\sim \mathcal{N}(0,{\sigma }^2)$，可得到第 $i$ 个目标变量 $t_i$ 的似然函数：

$p(t_i|\mathbf{w},{\mathbf{x}}_i)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2}{2{\sigma }^2}\right\}$

\qquad

• 一般假设观测误差 ${\epsilon }_i$ 满足独立同分布，可得到关于所有观测数据集 $\mathbf{X}$ 的目标值向量 $\mathbf{t}=[t_1,\cdots ,t_N{]}^T$ 的似然函数（参数为权向量 $\mathbf{w}$，$\mathbf{X}$ 表示数据集）：

$\begin{array}{rl}p(\mathbf{t}|\mathbf{w},\mathbf{X})& =p(t_1,\cdots ,t_N|\mathbf{w},\mathbf{X})\\ & =\prod _{i=1}^{N}p(t_i|\mathbf{w},{\mathbf{x}}_i)\\ & =\frac{1}{(\sqrt{2\pi }\sigma {)}^N}\exp \left\{-\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2\right\}\end{array}$

显然，此处的似然都是指 $p(\mathbf{t}|\mathbf{w},\mathbf{X})$ 随着参数 $\mathbf{w}$ 变化时的，而数据集 $\mathbf{X}$ 仅表示条件。

\qquad

2. 最大似然估计

\qquad 在得到了所有观测变量的联合概率密度 $p(\mathbf{t}|\mathbf{w},\mathbf{X})$ 之后，最直接的方法就是采用最大似然估计求取权向量 $\mathbf{w}$ 的解。

\qquad 对联合概率密度取对数，构成对数似然函数：

$\begin{array}{rl}\ln p(\mathbf{t}|\mathbf{w},\mathbf{X})& =\ln \frac{1}{(\sqrt{2\pi }\sigma {)}^N}-\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2\\ & =-\frac{N}{2}\ln (2\pi )-N\ln \sigma -\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2\end{array}$

观测数据集 $\mathbf{X}$ 出现在条件变量的位置，对于最大似然估计的求解而言可以忽略。

\qquad 采用最大似然估计求取作为参数的权向量 $\mathbf{w}$ 的值，也就是：

${\nabla }_{\mathbf{w}}\ln p(\mathbf{t}|\mathbf{w},\mathbf{X})=\frac{1}{{\sigma }^2}\sum _{i=1}^N\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]\mathbf{\varphi }({\mathbf{x}}_i)=0$

可以看出，在假设观测误差 $\epsilon$ 满足高斯噪声的情况下：
最大化对数似然函数 $\ln p(\mathbf{t}|\mathbf{w},\mathbf{X})$，实际上就是最小化平方和误差函数 $\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2$

$⟹$$\sum _{i=1}^N{t}_i\mathbf{\varphi }({\mathbf{x}}_i)=\sum _{i=1}^N[{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)]\mathbf{\varphi }({\mathbf{x}}_i)$

$⟹$$\underset{{\Phi }^T\mathbf{t}}{\underbrace{\left[\begin{array}{ccc}\mathbf{\varphi }({\mathbf{x}}_1)& \cdots & \mathbf{\varphi }({\mathbf{x}}_N)\end{array}\right]\left[\begin{array}{c}{t}_1\\ ⋮\\ t_N\end{array}\right]}}=\underset{{\Phi }^T\Phi \mathbf{w}}{\underbrace{\left[\begin{array}{ccc}\mathbf{\varphi }({\mathbf{x}}_1)& \cdots & \mathbf{\varphi }({\mathbf{x}}_N)\end{array}\right]\left[\begin{array}{c}\mathbf{\varphi }({\mathbf{x}}_1{)}^T\mathbf{w}\\ ⋮\\ \mathbf{\varphi }({\mathbf{x}}_N{)}^T\mathbf{w}\end{array}\right]}}$

$⟹$${\Phi }^T\mathbf{t}={\Phi }^T\Phi \mathbf{w}$

\qquad
\qquad 可求得最大似然解：

${\mathbf{w}}_{ML}=({\Phi }^T\Phi {)}^{-1}{\Phi }^T\mathbf{t}$

或将误差函数写成矩阵的形式 $(\mathbf{t}-\Phi \mathbf{w}{)}^T(\mathbf{t}-\Phi \mathbf{w})$，采用矩阵的微分得到

\qquad 其中，$\mathbf{\varphi }(\mathbf{x})=[{\varphi }_1(\mathbf{x}),\cdots ,{\varphi }_M(\mathbf{x}){]}^T$，$\mathbf{t}=[t_1,\cdots ,t_N{]}^T$

$\Phi =\left[\begin{array}{c}\mathbf{\varphi }({\mathbf{x}}_1{)}^T\\ \mathbf{\varphi }({\mathbf{x}}_2{)}^T\\ ⋮\\ \mathbf{\varphi }({\mathbf{x}}_N{)}^T\end{array}\right]=\left[\begin{array}{cccc}{\varphi }_1({\mathbf{x}}_1)& {\varphi }_2({\mathbf{x}}_1)& \cdots & {\varphi }_M({\mathbf{x}}_1)\\ {\varphi }_1({\mathbf{x}}_2)& {\varphi }_2({\mathbf{x}}_2)& \cdots & {\varphi }_M({\mathbf{x}}_2)\\ ⋮& ⋮& & ⋮\\ {\varphi }_1({\mathbf{x}}_N)& {\varphi }_2({\mathbf{x}}_N)& \cdots & {\varphi }_M({\mathbf{x}}_N)\end{array}\right]$

$\begin{array}{rl}{\Phi }^T\Phi & =\left[\begin{array}{cccc}\mathbf{\varphi }({\mathbf{x}}_1)& \mathbf{\varphi }({\mathbf{x}}_2)& \cdots & \mathbf{\varphi }({\mathbf{x}}_N)\end{array}\right]\left[\begin{array}{c}\mathbf{\varphi }({\mathbf{x}}_1{)}^T\\ \mathbf{\varphi }({\mathbf{x}}_2{)}^T\\ ⋮\\ \mathbf{\varphi }({\mathbf{x}}_N{)}^T\end{array}\right]\\ & =\sum _{i=1}^N\mathbf{\varphi }({\mathbf{x}}_i)\mathbf{\varphi }({\mathbf{x}}_i{)}^T\end{array}$

\qquad

3. 正则最小二乘(Regularized least-squares)估计

\qquad 上节在求最大似然解的时候，对数似然函数为：

$\ln p(\mathbf{t}|\mathbf{w},\mathbf{X})=-\frac{N}{2}\ln (2\pi )-N\ln \sigma -\frac{1}{{\sigma }^2}\left\{\frac{1}{2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2\right\}$
\qquad
\qquad 可以看出，对数似然函数的最后一项就是平方和误差函数：

$\begin{array}{rl}{E}_D(\mathbf{w})& =\frac{1}{2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2\\ & =\frac{1}{2}(\mathbf{t}-\Phi \mathbf{w}{)}^T(\mathbf{t}-\Phi \mathbf{w})\\ & =\frac{1}{2}\Vert \mathbf{t}-\Phi \mathbf{w}{\Vert }^2\end{array}$

\qquad 因此，最大化对数似然函数 $\ln p(\mathbf{t}|\mathbf{w},\mathbf{X})$，实际上就是最小化平方和误差函数 $E_D(\mathbf{w})$。为了防止以平方和误差函数作为代价函数时出现过拟合，可以添加正则化项来控制。因此，可以将代价函数修正为：

$E_D(\mathbf{w})+\lambda E_W(\mathbf{w})$

\qquad 其中，$\lambda$ 是正则化系数，用于平衡误差函数和正则化项的重要性。

\qquad 常用的正则化项可采用：

$E_W(\mathbf{w})=\frac{1}{2}\Vert \mathbf{w}{\Vert }^2=\frac{1}{2}{\mathbf{w}}^T\mathbf{w}$

\qquad 因此，得到了添加正则化项之后的代价函数 $F(\mathbf{w})$，也就是：

$F(\mathbf{w})=\frac{1}{2}\Vert \mathbf{t}-\Phi \mathbf{w}{\Vert }^2+\frac{\lambda }{2}\Vert \mathbf{w}{\Vert }^2$

$(1)$ 当 $\lambda =0$，表明对“由观测训练样本集 $\mathbf{X}$ 所描述的的观测模型”有完全的把握
$(2)$ 当 $\lambda =\infty$，表明对“由观测训练样本集 $\mathbf{X}$ 所描述的的观测模型”没有把握

\qquad
\qquad 为了最小化代价函数 $F(\mathbf{w})$，令 ${\nabla }_{\mathbf{w}}F(\mathbf{w})=-{\Phi }^T\mathbf{t}+\left({\Phi }^T\Phi +\lambda \mathbf{I}\right)\mathbf{w}=0$
　
\qquad 可得到正则最小二乘解：

$\mathbf{w}={\left[{\Phi }^T\Phi +\lambda \mathbf{I}\right]}^{-1}{\Phi }^T\mathbf{t}$

如果采用基本线性模型 $t=\sum _{j=1}^M{w}_j{x}_j+\epsilon$ ，满足 $\mathbf{\varphi }({\mathbf{x}}_i)={\mathbf{x}}_i$
$t_i=\sum _{j=1}^M{w}_j{x}_{ij}+{\epsilon }_i$　或者　$t_i={\mathbf{w}}^T{\mathbf{x}}_i+{\epsilon }_i$
平方和误差：$E_D(\mathbf{w})=\sum _{i=1}^N{\epsilon }_i^2=\sum _{i=1}^N(t_i-{\mathbf{w}}^T{\mathbf{x}}_i{)}^2=(\mathbf{t}-\Phi \mathbf{w}{)}^T(\mathbf{t}-\Phi \mathbf{w})$
　
此时，$\Phi =\left[\begin{array}{c}\mathbf{\varphi }({\mathbf{x}}_1{)}^T\\ \mathbf{\varphi }({\mathbf{x}}_2{)}^T\\ ⋮\\ \mathbf{\varphi }({\mathbf{x}}_N{)}^T\end{array}\right]=\left[\begin{array}{c}{\mathbf{x}}_1^T\\ {\mathbf{x}}_2^T\\ ⋮\\ {\mathbf{x}}_N^T\end{array}\right]=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1M}\\ x_{21}& x_{22}& \cdots & x_{2M}\\ ⋮& & & \\ x_{N1}& x_{N2}& \cdots & x_{NM}\end{array}\right]$
　
采用正则化方法防止过拟合，优化函数可表示为：$F(\mathbf{w})=(\mathbf{t}-\Phi \mathbf{w}{)}^T(\mathbf{t}-\Phi \mathbf{w})+\lambda {\mathbf{w}}^T\mathbf{w}$
　
令 ${\nabla }_{\mathbf{w}}F(\mathbf{w})=0$，可得到：$\mathbf{w}=[{\Phi }^T\Phi +\lambda \mathbf{I}{]}^{-1}{\Phi }^T\mathbf{t}$

\qquad

4. 最大后验估计

\qquad 从贝叶斯分析的观点，需要考虑作为参数的权向量 $\mathbf{w}=[w_1,\cdots ,w_M{]}^T$ 的先验概率。
\qquad

• 仍然假设权向量的各元素 $w_i$ 满足独立同分布，且满足 $w_i\sim \mathcal{N}(0,{\sigma }_w^2)$，可得到先验概率：

$\begin{array}{rl}p(\mathbf{w})& =p(w_1,\cdots ,w_M)\\ & =\prod _{i=1}^{M}p(w_i)\\ & =\frac{1}{(\sqrt{2\pi }{\sigma }_w{)}^M}\prod _{i=1}^M\exp \left(-\frac{w_i^2}{2{\sigma }_w^2}\right)\\ & =\frac{1}{(\sqrt{2\pi }{\sigma }_w{)}^M}\exp \left(-\frac{1}{2{\sigma }_w^2}\sum _{i=1}^M{w}_i^2\right)\\ & =\frac{1}{(\sqrt{2\pi }{\sigma }_w{)}^M}\exp \left(-\frac{\Vert \mathbf{w}{\Vert }^2}{2{\sigma }_w^2}\right)\end{array}$
\qquad

• 由贝叶斯公式，可得到后验概率满足：

$p(\mathbf{w}|\mathbf{t},\mathbf{X})=\frac{p(\mathbf{t}|\mathbf{w},\mathbf{X})p(\mathbf{w})}{p(\mathbf{t})}\propto p(\mathbf{t}|\mathbf{w},\mathbf{X})p(\mathbf{w})$
\qquad
\qquad 采用第 $1$ 部分描述的线性回归的概率模型，目标向量 $\mathbf{t}$ 的似然函数为：

$p(\mathbf{t}|\mathbf{w},\mathbf{X})=\frac{1}{(\sqrt{2\pi }\sigma {)}^N}\exp \left\{-\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2\right\}$

\qquad 那么

$\begin{array}{rl}p(\mathbf{w}|\mathbf{t},\mathbf{X})& \propto \frac{1}{(\sqrt{2\pi }\sigma {)}^N}\exp \left\{-\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2\right\}\frac{1}{(\sqrt{2\pi }{\sigma }_w{)}^M}\exp \left(-\frac{\Vert \mathbf{w}{\Vert }^2}{2{\sigma }_w^2}\right)\\ & \propto \exp \left\{-\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2-\frac{\Vert \mathbf{w}{\Vert }^2}{2{\sigma }_w^2}\right\}\\ & \propto -\left\{\frac{1}{2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2+\frac{\lambda }{2}\Vert \mathbf{w}{\Vert }^2\right\}\end{array}$

\qquad\qquad 其中，定义 $\lambda =\frac{{\sigma }^2}{{\sigma }_w^2}$
\qquad

• 权向量 $\mathbf{w}$ 的最大后验估计值为：

$\begin{array}{rl}{\mathbf{w}}_{MAP}& =\underset{\mathbf{w}}{\max }p(\mathbf{w}|\mathbf{t},\mathbf{X})\\ & =\underset{\mathbf{w}}{\max }\left\{-\frac{1}{2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2-\frac{\lambda }{2}\Vert \mathbf{w}{\Vert }^2\right\}\end{array}$
\qquad
\qquad 相当于定义了代价函数 $F(\mathbf{w})$，也就是：

$\begin{array}{rl}F(\mathbf{w})& =\frac{1}{2}\sum _{i=1}^N{\left[t_i-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_i)\right]}^2+\frac{\lambda }{2}\Vert \mathbf{w}{\Vert }^2\\ & =\frac{1}{2}\Vert \mathbf{t}-\Phi \mathbf{w}{\Vert }^2+\frac{\lambda }{2}\Vert \mathbf{w}{\Vert }^2\end{array}$

\qquad
\qquad 最大后验估计值就为： ${\mathbf{w}}_{MAP}=\underset{\mathbf{w}}{\min }F(\mathbf{w})$

\qquad 显然，在假设权向量各元素满足独立同分布、且 $w_i\sim \mathcal{N}(0,{\sigma }_w^2)$ 的条件下，权向量 $\mathbf{w}$ 的最大后验解，就是第 $3$ 节所描述的正则最小二乘解。
\qquad
\qquad 此时，权向量 $\mathbf{w}$ 的最大后验解为：$\mathbf{w}={\left[{\Phi }^T\Phi +\lambda \mathbf{I}\right]}^{-1}{\Phi }^T\mathbf{t}$

\qquad 特别地，当取 $\lambda =0({\sigma }_w^2=\infty )$ 时，权向量各元素 $w_i$ 近似于均匀分布，也就是完全忽视权向量 $\mathbf{w}$ 的先验信息，此时的最大后验解就等价于最大似然解，也就是：${\mathbf{w}}_{MAP}={\mathbf{w}}_{ML}={\left({\Phi }^T\Phi \right)}^{-1}{\Phi }^T\mathbf{t}$。
\qquad

举例：双月数据分类

import numpy as np
import matplotlib.pyplot as plt

def gen_lineardata(weight,interval):
    y = -(weight[0]*interval + weight[2])/weight[1]
    return y

def halfmoon(rad, width, dist, n_samp):      
    if n_samp%2 != 0:  
        n_samp += 1      
    data = np.zeros((3,n_samp))      
    rd = np.random.random((2,n_samp//2))  
    radius = (rad-width//2) + width*rd[0,:] 
    theta = np.pi*rd[1,:]          
    x1     = radius*np.cos(theta)  
    y1     = radius*np.sin(theta) + dist/2  
    label1 = np.ones((1,len(x1)))           # label= 1 for Class 1  
    rd = np.random.random((2,n_samp//2))   
    radius = (rad-width//2) + width*rd[0,:] 
    theta = np.pi*rd[1,:]
    x2    = radius*np.cos(-theta) + rad  
    y2    = radius*np.sin(-theta) - dist/2  
    label2= -1*np.ones((1,len(x2)))           # label= 0 for Class 2       
    data[0,:]=np.concatenate([x1,x2])
    data[1,:]=np.concatenate([y1,y2])
    data[2,:]=np.concatenate([label1,label2],axis=1)    
    shuffle_seq = np.random.permutation(np.arange(n_samp))  
    data_shuffle = data[:,shuffle_seq]
    return data,data_shuffle

def RLS(xhat,target,lambda0):
    Phi = np.asmatrix(xhat)
    t = np.asmatrix(target)
    print(Phi.T*Phi)
    print(Phi.T*Phi+lambda0*np.eye(len(xhat.T)))    
    return np.array((Phi.T*Phi+lambda0*np.eye(len(xhat.T))).I*Phi.T*t)    

if __name__ == "__main__":
    dNum = 800    
    data,data_shuffle = halfmoon(10,6,1,dNum)
    #data,data_shuffle = halfmoon(10,6,-4,dNum)
    pos_data = data[:,0:dNum//2]
    neg_data = data[:,dNum//2:dNum]  
    training_data = data_shuffle.T
    tmp1 = training_data[0:dNum,0:2]
    tmp2 = np.ones((dNum,1))
    xhat = np.concatenate((tmp1,tmp2),axis=1)
    target = training_data[0:dNum,2:]    
    interval = np.linspace(-12,20,100)    
    weight = RLS(xhat,target,0)
    print('RLS:',weight.flatten())    
    y = gen_lineardata(weight,interval)
    plt.figure   
    plt.plot(interval,y,'k')       
    plt.plot(pos_data[0,:],pos_data[1,:],'b+')
    plt.plot(neg_data[0,:],neg_data[1,:],'r+')        
    plt.title('Regularized least squares')     
    plt.show()

\qquad

距离=1，半径=10，宽度=6

\qquad

距离=-4，半径=10，宽度=6