Typesetting math: 100%

PRML-3.3 贝叶斯线性回归

记号说明

1.X={x1,...,xN}N,{xn},n=1,2,...,N,x_train,D
2.t={t1,...,tn},y_train
3.t=y(x),xt
4.p(t|x)
5.线y(x,w)=w0+w1x1+...+wDxD=w0+j=1M1wjϕ(x)=wTϕ(x)
6.M
7.ϕj(x)(basis function)j=0,...,M1
8.w0 bias parameter
9.w=(w0,...,wM1)T
10.ϕ(x)=(ϕ0,...,ϕM1)T
11.ϕj(x)=xj,
12.β,,p(t|x,w,β)=N(t|y(x,w),β1)
13.Φ=(ϕ0(x1)ϕ1(x1)ϕM1(x1)ϕ0(x2)ϕ1(x2)ϕM1(x2)ϕ0(xN)ϕ1(xN)ϕM1(xN))N×M
14.Φ+(ΦTΦ)1ΦT,摩尔彭罗斯伪逆(Moore-Penrose pseudo-inverse)
15.ED(w)=12n=1N{tnw0j=1M1wjϕj(xn)}2=12n=1N{tnwTϕ(xn)}2,
15.EW(w)=12wTw,

3.3.1

引入模型参数w的先验概率分布,吧噪声精度参数β当做已知参数
于是对应的共轭先验是形式为

(3.48)p(w)=N(w|m0,S0)

均值和方差分别为$ m_0, S_0 (S_0$是已知参数)
由于共轭的性质,后验也是高斯分布

(3.49)p(w|t)=N(w|mN,SN)

其中

(3.50)mN=SN(S01m0+βΦTt)(3.51)SN1=S01+βΦTΦ


证明
贝叶斯定理p(w|t)p(t|w)p(w)
代入3.10 ,3.48
p(t|X,w,β)=n=1NN(tn|wTϕ(xn),β1)
p(w)=N(w|m0,S0)

p(w|t)exp[β2(t1wTϕ1)2]exp[β2(t2wTϕ2)2]...exp[β2(tnwTϕn)2]exp[12(wm0)TS01(wm0)]
=exp[β2{(t1wTϕ1)2+...+(tnwTϕn)2}]exp[12(wm0)TS01(wm0)]
=exp[β2(t1wTϕ1...tnwTϕn)T(t1wTϕ1...tnwTϕn)]exp[12(wm0)TS01(wm0)]
=exp[β2(TΦw)T(TΦw)]exp[12(wm0)TS01(wm0)]
其中
T=(t1...tn)TΦwT=(ϕ1Tw...ϕnTw)=(wTϕ1...wTϕn)T
=exp[12(βTTTβTTΦwβ(Φw)TT+βwTΦTΦTw+wTS01wm0TS01wwTS01m0+m0TS01m0)]
=exp[12(wT(S01+βΦTΦ)wwT(S01m0+βΦTT)(βTTΦ+m0TS01)w+m0TS01m0+βTTT)],同色的相结合
3.50式mN=SN(S01m0+βΦTT)
3.51式SN1=S01+βΦTΦ
用配方法,不确定可以反向乘开来验证
=exp[12(wmN)TSN1(wmN)]exp[12(m0TS01m0+βTT+mNTSNTmN)]
=C×exp[12(wmN)TSN1(wmN)],因为后面一项和w没有关系


方法2 利用前面2.3.3 章节的结论

对于$ x y x $的条件高斯分布:

(2.113)p(x)=N(x|μ,Λ1)

(2.114)p(y|x)=N(y|Ax+b,L1)

那么$ y x y $的条件高斯分布为:

(2.115)p(y)=N(y|Aμ+b,L1+AΛ1AT)

(2.116)p(x|y)=N(x|Σ{ATL(yb)+Λμ},Σ)

其中

(2.117)Σ=(Λ+ATLA)1


0均值的先验

(3.53)mN=βSNΦTt(3.54)SN1=αI+βΦTΦ


后验分布的对数是由对数似然与先验的对数的和给出的一个关于w的函数,形式为:

(3.55)lnp(w|t)=β2n=1N{tnwTϕ(xn)}2α2wTw+const

关于w来最大化这个后验分布等价于最小化加上一个二次正则项的平方和误差函数,正则项对应λ=α/β的式(3.27)。


贝叶斯线性回归的顺序学习

看图说话
0.模型选用简单的y(x,w)=w0+w1x
1.如果数据点时顺序到达的,那么任何一个阶段的后验概率分布都可以看成后续数据点的先验
2.假设先验的β和后验的α都是已知的
3.当新数据点被观测到的时候,当前的后验分布变成了先验分布
4.第一行对应于未观测到任何数据点之前的情况,因为w的先验假设是正态的,所以第一行中间一列的图像是正态的,第一行右边一列是6个y_test的值,用公式y(x,w)=w0+w1x,其中w0,w1都是先验采样出来的,也就是这一格都是正态分布采样出来的,看起来杂乱无章
5.左中两列,都有一个白色的+号,表示的是w0,w1的真实数据
6.左列是似然函数,一行二列的先验 × 二行一列的似然 = 二行二列的后验,下面的都是以此类推
7.第二行是观测到第一个观测数据后的更新后验的情况
8.第三行是观测到第二个观测后的情况
9.观测到20个数据的情况


证明-习题3.8

Solution.
ΦN=[ϕ1T;;ϕNT],tN=[t1,,tN]T,则已知 N 个样本的后验估计可以表示为

mN=SN(S01m0+βΦNTtN)SN1=S01+βΦNTΦN

ΦN+1=[ΦN;ϕN+1],tN+1=[tNT,tN+1]T,则

ΦN+1TΦN+1=[ΦNT,ϕN+1][ΦNϕN+1]=ΦNTΦNT+ϕN+1ϕN+1TΦN+1TtN+1=[ΦNT,ϕN+1][tNtN+1]=ΦNTtN+tN+1ϕN+1

SN+11=S01+βΦN+1TΦN+1=S01+βΦNTΦN+βϕN+1ϕN+1T=SN1+βϕN+1ϕN+1T

另外

mN+1=SN+1(S01m0+βΦN+1TtN+1)=SN+1(S01m0+βΦNTtN+βtN+1ϕN+1)=SN+1(SN1mN+βtN+1ϕN+1)

即得到了增加第 N+1 个样本时的更新公式。可以看到,递推式和通项公式形式上十分接近。


其他形式的参数先验也可以被考虑。例如,我们可以推广高斯先验分布,得到:

(3.56)p(wα)=[q2(α2)1q1Γ(1q)]Mexp(α2j=0M1|wj|q)

q=2时对应着高斯分布,并且只有在这种情形下的先验才与式(3.10)给出的似然函数共轭。找到w的后验概率分布的最大值等价于找到正则化误差函数(3.29)的最小值。在高斯先验的情况下,后验概率分布的众数等于均值,但是当q2时这个性质就不成立了。


python实现

点击查看代码


import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt
%matplotlib inline

from prml.preprocess import GaussianFeature, PolynomialFeature, SigmoidalFeature
from prml.linear import (
    BayesianRegression,
    EmpiricalBayesRegression,
    LinearRegression,
    RidgeRegression
)

np.random.seed(1234)

def create_toy_data(func, sample_size, std, domain=[0, 1]):
    x = np.linspace(domain[0], domain[1], sample_size)
    np.random.shuffle(x)
    t = func(x) + np.random.normal(scale=std, size=x.shape)
    return x, t

# 自定义图3.7
def linear(x):
    return -0.3 + 0.5 * x


x_train, y_train = create_toy_data(linear, 20, 0.1, [-1, 1])
x = np.linspace(-1, 1, 100)
w0, w1 = np.meshgrid(
    np.linspace(-1, 1, 100),
    np.linspace(-1, 1, 100))
w = np.array([w0, w1]).transpose(1, 2, 0)

feature = PolynomialFeature(degree=1)
X_train = feature.transform(x_train)
X = feature.transform(x)
model = BayesianRegression(alpha=1., beta=100.)

for begin, end in [[0, 0], [0, 1], [1, 2], [2, 3], [3, 20]]:
    model.fit(X_train[begin: end], y_train[begin: end])
    plt.subplot(1, 2, 1)
    plt.scatter(-0.3, 0.5, s=200, marker="x")
    plt.contour(w0, w1, multivariate_normal.pdf(w, mean=model.w_mean, cov=model.w_cov))
    plt.gca().set_aspect('equal')
    plt.xlabel("$w_0$")
    plt.ylabel("$w_1$")
    plt.title("prior/posterior")

    plt.subplot(1, 2, 2)
    plt.scatter(x_train[:end], y_train[:end], s=100, facecolor="none", edgecolor="steelblue", lw=1)
    plt.plot(x, model.predict(X, sample_size=6), c="orange")
    plt.xlim(-1, 1)
    plt.ylim(-1, 1)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.show()

点击查看代码

import numpy as np
from prml.linear.regression import Regression


class BayesianRegression(Regression):
    """
    贝叶斯回归模型

    w ~ N(w|0, alpha^(-1)I)
    y = X @ w
    t ~ N(t|X @ w, beta^(-1))
    """

    def __init__(self, alpha: float = 1., beta: float = 1.):
        self.alpha = alpha
        self.beta = beta
        self.w_mean = None
        self.w_precision = None

    def _is_prior_defined(self) -> bool:
        return self.w_mean is not None and self.w_precision is not None

    def _get_prior(self, ndim: int) -> tuple:  # -> tuple 表示返回一个tuple
        '''
        :param ndim:
        :return:  w_mean  先验的均值
                    w_precision 先验的精度
        '''
        if self._is_prior_defined():
            return self.w_mean, self.w_precision
        else:
            return np.zeros(ndim), self.alpha * np.eye(ndim)

    def fit(self, X: np.ndarray, t: np.ndarray):
        """
        给定训练集,更新贝叶斯参数

        Parameters
        ----------
        X : (N, n_features) np.ndarray
            训练集输入
        t : (N,) np.ndarray
            训练集输出
        """

        # 先验的均值,先验的精度
        mean_prev, precision_prev = self._get_prior(
            np.size(X, 1))  # mean_prev =m_0 precision_prev =S_0= \alpha I 见PRML 3.48

        w_precision = precision_prev + self.beta * X.T @ X  # 这里是计算 w_precision=S_N^{-1}=S_0^{-1}+\beta X^T X,PRML 3.54

        # 计算m(x) 也就是后验的均值
        w_mean = np.linalg.solve(
            w_precision,
            precision_prev @ mean_prev + self.beta * X.T @ t
            # w_mean=m_N ,解方程 S_N^{-1} w_mean =S_0^{-1} * m_0 +\beta X^T t ,PRML 3.53
        )
        self.w_mean = w_mean
        self.w_precision = w_precision
        self.w_cov = np.linalg.inv(self.w_precision)  # w_cov = S_N , 协方差矩阵

    def predict(self, X: np.ndarray, return_std: bool = False, sample_size: int = None):
        """
        return mean (and standard deviation) of predictive distribution

        Parameters
        ----------
        X : (N, n_features) np.ndarray
            independent variable
        return_std : bool, optional
            flag to return standard deviation (the default is False)
        sample_size : int, optional
            number of samples to draw from the predictive distribution
            (the default is None, no sampling from the distribution)

        Returns
        -------
        y : (N,) np.ndarray
            mean of the predictive distribution
        y_std : (N,) np.ndarray
            standard deviation of the predictive distribution
        y_sample : (N, sample_size) np.ndarray
            samples from the predictive distribution
        """

        # 得到拟合的正态分布后,对这个正态分布采样
        if sample_size is not None:
            w_sample = np.random.multivariate_normal(
                self.w_mean, self.w_cov, size=sample_size
            )
            y_sample = X @ w_sample.T
            return y_sample
        y = X @ self.w_mean
        if return_std:
            y_var = 1 / self.beta + np.sum(X @ self.w_cov * X, axis=1)
            y_std = np.sqrt(y_var)
            return y, y_std
        return y


3.3.2 预测分布

预测分布如下

公式3.57


(3.57)p(t|T,α,β)=p(t|w,β)p(w|T,α,β)dw

1.p(t,T,α,β)=p(t,T,α,β,w)dw
=p(t|T,α,β,w)p(T,α,β,w)dw
=p(t|T,α,β,w)p(w|T,α,β)p(T,α,β)dw

2.p(t,T,α,β)=p(t|T,α,β)p(T,α,β)
1=2p(x,X,T)w
p(t|T,α,β)=p(t|T,α,β,w)p(w|T,α,β)dw

p(t|T,α,β,w)tα,T,p(t|T,α,β,w)=p(t|β,w)

p(t|T,α,β)=p(t|w,β)p(w|T,α,β)dw


目标变量的条件分布由式(3.8)给出,后验分布由式(3.49)给出。

目标变量的条件分布:

(3.8)p(t|x,w,β)=N(t|y(x,w),β1)

后验分布
(3.49)p(w|t)=N(w|mN,SN)

得到的预测分布的形式为:

(3.58)p(t|x,t,α,β)=N(t|mNTϕ(x),σN2(x))

其中预测分布的方差σN2(x)

(3.59)σN2(x)=1β+ϕ(x)TSNϕ(x)


当额外的数据点被观测到的时候,后验概率分布会变窄,从而可以证明σN+12(x)σN2(x)
证明-习题3.11

Hint.

ΦN=[ϕ0T;;ϕNT]

SN+11=S01+βΦN+1TΦN+1=S01+β[ΦNT,ϕN+1][ΦNϕN+1T]=S01+β(ΦNTΦN+ϕN+1ϕN+1T)=SN1+vvT(v=βϕN+1)

根据提示(3.110恒等式)(M+vvT)1=M1(M1v)(vTM1)1+vTM1v(M+vvT)1M1,令M=SN1,则有SN+1SN,则有σN+12σN2

3.3.3 贝叶斯模型的等价核解释

预测均值E(y_test)与训练集目标值y_train之间的线性关系3.60 3.61

y(x,mN)=n=1Nβϕ(x)TSNϕ(xn)tn=n=1Nk(x,xn)tn

其中函数
(3.62)k(x,x)=βϕ(x)TSNϕ(x)
被称为平滑矩阵(smoother matrix)或等价核(equivalent kernel)。像这样,通过训练集的目标值的线性组合来进行预测的回归函数被称为线性平滑(linear smoother)

用核函数表示线性回归是解决回归问题的另一种方法。在给定观测数据集的条件下,直接定义一个局部的核函数,而不是引入一组基函数来隐式地定义了一个等价的核,然后使用这个核函数对新的输入变量$ x 6.4Gaussianprocess x $值做预测。可以证明权值的和为1,即:

(3.64)n=1Nk(x,xn)=1

对于所有x成立。

posted @   筷点雪糕侠  阅读(252)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· spring官宣接入deepseek,真的太香了~
点击右上角即可分享
微信分享提示