[机器学习] 使用极大似然估计与贝叶斯估计拟合线性基函数

一、概率模型

概率模型可以写作：
$$P(\cdot |\theta )$$
指模型的分布类型已知（二项分布、正态/高斯分布，Bate分布等），但是模型的参数未知。基于独立同分布的样本数据D对模型参数进行估计。常用的方法有极大似然估计（MLE）、极大后验估计（MAP）、贝叶斯估计。

二、极大似然估计Maximum Likelihood Estimation

MLE的实质为计算在参数下D出现的概率，需要计算D中每个样本出现概率的乘积，即似然函数。表达式如下：
D = { x i } i = 1 n L i k e l i h o o d ( θ ) = P ( D ∣ θ ) = Π i = 1 n P ( x i ∣ θ ) 似然函数 l n L ( θ ) = Σ i = 1 n l n P ( x i ∣ θ ) 对数似然函数 D = \{x_i\}_{i=1}^n \\ Likelihood(\theta) = P(D|\theta) = \Pi_{i = 1}^nP(x_i|\theta) \quad 似然函数 \\ lnL(\theta) = \Sigma_{i=1}^nln P(x_i|\theta)\quad 对数似然函数 D={xi​}i=1n​Likelihood(θ)=P(D∣θ)=Πi=1n​P(xi​∣θ)似然函数lnL(θ)=Σi=1n​lnP(xi​∣θ)对数似然函数
对似然函数取对数可以将概率乘法变为加法，将参数估计问题变为函数优化问题，从而可以得到能够使观察数据D出现概率最大的概率模型P参数来。

三、贝叶斯估计Bayesian

贝叶斯估计表达式如下。其中p(A|B)为B条件下A的概率，w为参数，D为独立同分布数据样本：
$$p(w|D)=\frac{p(D|w)p(w)}{p(D)}$$
贝叶斯模型生成参数模型的概率分布（后验概率），其实质为基于样本数据调整先验概率，得到更合理的后验概率。

四、线性基函数模型

用于拟合非线性数据，普通模型的计算为：
$$y=x^t\cdot w$$
线性基函数即对输入x进行线性变换后，再与参数进行计算。一般是为了对输入特征的维度进行变换，参数维度K决定了模型的描述能力：K=0，模型为一条水平线；K=1，模型为一条直线…然而K的值不是越高越好，过大的K会使得模型容易过拟合。

$$\begin{gathered} y={\varphi }^t(x)\cdot w+ϵ(x\in R^D,w\in R^K,ϵ\sim N(0,{\sigma }^2)) \\ \varphi (x)=\left|\begin{array}{c}{\varphi }_0(x)\\ {\varphi }_1(x)\\ ⋮\\ {\varphi }_{K-1}(x)\end{array}\right|=\left|\begin{array}{c}1\\ x\\ x^2\\ ⋮\\ x^{K-1}\end{array}\right|\in R^K \\ ⟹y=w_{0}1+w_{1}x+\cdots +w_{k-1}{x}^{k-1}+ϵ \end{gathered}$$

五、实验

1、训练集生成

• 使用sin函数和方差为0.2高斯噪声生成S、M、L三个数据集，大小分别为10、30、60。

# 使用sin函数和方差为2高斯噪声生成非线性数据样本
import numpy as np
import random

# n 生成样本数目
def data_gen(num):
    # 在[-5,5]内生成num个样本
    X = np.linspace(-5,5,num,endpoint=True)
    Y = np.sin(X)

    # 加入高斯噪声
    for i in range(X.size):
        X[i] += random.gauss(0,0.2)
        Y[i] += random.gauss(0,0.2)
    return X,Y

2、模型训练

• 使用线性基函数作为回归模型，针对不同的K，使用MLE在不同数据集上进行拟合。

对于观测数据X和Y，有：
$$\begin{gathered} w_{MLE}=({\Phi }^T\Phi {)}^{-1}{\Phi }^{T}y \\ {\sigma }_{MLE}^2=\frac{1}{N}\sum _{i=1}^N(y_i-{\Phi }^T(x_i)w{)}^2 \end{gathered}$$
其中：
$$\Phi =\left|\begin{array}{ccc}{\varphi }_0(x_1)& \cdots & {\varphi }_{K-1}(x_1)\\ {\varphi }_0(x_2)& \cdots & {\varphi }_{K-1}(x_2)\\ ⋮& \ddots & ⋮\\ {\varphi }_0(x_N)& \cdots & {\varphi }_{K-1}(x_N)\end{array}\right|\in R^{N\ast K}$$

import math

# 线性基函数模型，使用MLE拟合
class model_MLE:
    # K为w的维数，w为一个K维向量，sigma为噪声方差，XY为训练数据
    def __init__(self,K,X,Y):
        self.K = K
        
        self.X = X
        self.Y = Y
        
        self.Phi = np.array(self.Phi_gen())
        self.w = self.w_gen()
        self.sigma = self.sigma_gen()
    
    # 极大似然估计计算w
    def w_gen(self):
        w = np.dot(np.dot(np.linalg.inv(np.dot(self.Phi.T,self.Phi)),self.Phi.T),self.Y)
        return w
    
    # 极大似然估计计算sigma
    def sigma_gen(self):
        sigma = 0
        for i in range(self.X.size):
            sigma += math.pow((self.Y[i] - np.dot(np.array(self.Phi_func(self.X[i])).T, self.w)),2)
        sigma /= self.X.size
        return sigma
    
    # 依据输入X和K生成Phi矩阵
    def Phi_gen(self):
        Phi = []
        for i in range(self.X.size):
            Phi_x = self.Phi_func(self.X[i])
            Phi.append(Phi_x)
        return Phi
    
    # Phi函数，对单个x进行处理
    def Phi_func(self,x):
        # 对x进行维度变换
        Phi_x =[]
        item = 1
        for i in range(self.K):
            Phi_x.append(item)
            item *= x
#         Phi_x = np.array(Phi_x)
        return Phi_x

    def forcast(self,x):
        Phi_x = self.Phi_func(x)
#         print(Phi_x,self.w)
        y = np.dot(Phi_x,self.w)
        return y

• 使用线性基函数作为回归模型，针对不同的K，使用贝叶斯估计在不同数据集上进行拟合。

贝叶斯估计将模型参数w和噪声视为一个参数未知的概率分布，通过求该概率分布的均值和方差就可以得到模型w的分布概率。设模型的参数的先验分布概率服从高斯分布：（其中I为单位矩阵）
$$\begin{gathered} P(w)\sim N({\mu }_0,{\Sigma }_0)=N(0,\frac{1}{\lambda }I) \\ ϵ\sim N(0,{\sigma }^2) \end{gathered}$$
有：
$$\begin{gathered} P(w|y,X,{\sigma }^2)\sim N({\mu }_{new}|{\Sigma }_{new}) \\ {\Sigma }_{new}=(\frac{1}{{\sigma }^2}{\Phi }^T\Phi +{\Sigma }_0^{-1}{)}^{-1}=(\frac{1}{{\sigma }^2}{\Phi }^T\Phi +\lambda I{)}^{-1} \\ {\mu }_{new}={\Sigma }_{new}(\frac{1}{{\sigma }^2}{\Phi }^{T}y+{\Sigma }_0^{-1}{\mu }_0)={\Sigma }_{new}(\frac{1}{{\sigma }^2}{\Phi }^{T}y) \end{gathered}$$

# 线性基函数模型参数概率分布，使用贝叶斯估计拟合，得到模型参数的后验概率分布
class distribution_BYS:
    # K为w的维数，w为一个K维向量，sigma为噪声的方差
    def __init__(self,K,X,Y):
        self.K = K
        
        self.X = X
        self.Y = Y
        
        # 模型参数先验分布
        self.mu = 0
#       I为单位矩阵
        self.I = np.identity(self.K, dtype = int)
        self.sigma = 0.25*self.I
        
        # 模型参数后验分布
        self.Phi = self.Phi_gen()
        self.sigma_new = self.sigma_new_gen()
        self.mu_new = self.mu_new_gen()
        
    def sigma_new_gen(self):
        sigma_new = np.linalg.inv(np.dot(self.Phi.T,self.Phi)/0.2+4*self.I)
        return sigma_new
    
    def mu_new_gen(self):
        mu_new = np.dot(self.sigma_new,(np.dot(self.Phi.T,self.Y)/0.2))
        return mu_new
    
     # 依据输入X和K生成Phi矩阵
    def Phi_gen(self):
        Phi = []
        for i in range(self.X.size):
            Phi_x = self.Phi_func(self.X[i])
            Phi.append(Phi_x)
        Phi = np.array(Phi)
        return Phi
    
    # Phi函数，对单个x进行处理
    def Phi_func(self,x):
        # 对x进行维度变换
        Phi_x =[]
        item = 1
        for i in range(self.K):
            Phi_x.append(item)
            item *= x
#         Phi_x = np.array(Phi_x)
        return Phi_x
    
    # 基于模型参数后验概率分布随机生成w，实现模型的随机取样
    def model_gen(self):
        w = []
        for i in range(self.K):
            w.append(random.gauss(self.mu_new[i],self.sigma_new[i][i]))
        w = np.array(w)
        return w

• 类实例化，使用不同的K在不同大小的数据集上进行拟合，代码略。

3、可视化

• 样本数据可视化：

# 样本数据可视化
import matplotlib.pyplot as plt

plt.figure(figsize = (12,8))

plt.subplot(2,2,1)
plt.plot(S_X,S_Y,linestyle='',marker='.',label = 'S')
plt.plot(np.linspace(-5,5),np.sin(np.linspace(-5,5)),label = 'Sin(x)')
plt.legend(loc = 1)

plt.subplot(2,2,2)
plt.plot(M_X,M_Y,linestyle='',marker='.',label = 'M')
plt.plot(np.linspace(-5,5),np.sin(np.linspace(-5,5)),label = 'Sin(x)')
plt.legend(loc = 1)

plt.subplot(2,2,3)
plt.plot(L_X,L_Y,linestyle='',marker='.',label = 'L')
plt.plot(np.linspace(-5,5),np.sin(np.linspace(-5,5)),label = 'Sin(x)')
plt.legend(loc = 1)

plt.show()

• MLE拟合结果可视化：

# print(model_MLE_1.forcast(X[0]))
# MLE拟合结果可视化
def model_show(model,lab):
    x = np.linspace(-5,5,100)
    y = []
    for i in range(x.size):
        y.append(model.forcast(x[i]))
    plt.plot(x,y,label = lab)

plt.figure(figsize = (12,8))

plt.subplot(2,2,1)
plt.plot(S_X,S_Y,linestyle='',marker='.',label = 'S')
model_show(model_MLE_1,"K=2")
model_show(model_MLE_2,"K=4")
model_show(model_MLE_3,"K=8")
plt.legend(loc = 1)

plt.subplot(2,2,2)
plt.plot(M_X,M_Y,linestyle='',marker='.',label = 'M')
model_show(model_MLE_4,"K=2")
model_show(model_MLE_5,"K=4")
model_show(model_MLE_6,"K=8")
plt.legend(loc = 1)

plt.subplot(2,2,3)
plt.plot(L_X,L_Y,linestyle='',marker='.',label = 'L')
model_show(model_MLE_7,"K=2")
model_show(model_MLE_8,"K=4")
model_show(model_MLE_9,"K=8")
plt.legend(loc = 1)

plt.show()

• 贝叶斯估计拟合结果可视化：

# 贝叶斯拟合结果可视化
def dis_show(dis):
    x = np.linspace(-5,5,100)
    for i in range(20):
        w = dis.model_gen()
#         print(w)
        y = []
        for j in range(x.size):
            y.append(np.dot(dis.Phi_func(x[j]),w.T))
        plt.plot(x,y)


def draw(loc,dis,data_lab,tit):
    plt.subplot(3,3,loc)
    if data_lab == 'S':
        plt.plot(S_X,S_Y,linestyle='',marker='.',label = data_lab)
    elif data_lab == 'M':
        plt.plot(M_X,M_Y,linestyle='',marker='.',label = data_lab)
    else:
        plt.plot(L_X,L_Y,linestyle='',marker='.',label = data_lab)
    dis_show(dis)
    plt.title(tit)
    plt.legend(loc = 1)

plt.figure(figsize = (12,9))
draw(1,dis_BYS_1,'S','K=2')
draw(2,dis_BYS_2,'S','K=4')
draw(3,dis_BYS_3,'S','K=8')

draw(4,dis_BYS_4,'M','K=2')
draw(5,dis_BYS_5,'M','K=4')
draw(6,dis_BYS_6,'M','K=8')

draw(7,dis_BYS_7,'L','K=2')
draw(8,dis_BYS_8,'L','K=4')
draw(9,dis_BYS_9,'L','K=8')

plt.show()