3.1线 性 回 归
线性回归通过拟合关于观测数据的线性方程来建模两个变量间的关系，其中一个变量是自变量，另一个是因变量，自变量与因变量都可以是多元变量。多元自变量对应事物的多个输入特征，多元因变量对应事物的多种输出信息。例如，可以通过一个线性回归模型来关联熊猫的性别、年龄、体重和食量,其中性别、年龄和体重作为自变量，食量作为因变量。
下面给出线性回归的原理与表示。给定有N个样本的数据集 ={yi,xi1,xi2,…,xiD}Ni=1，线性回归模型假设因变量yi与自变量xi（由{xi1,xi2,…,xiD}构成的D维向量）间是线性关系。此关系通过回归系数β构建,模型的形式为:yi=β01+β1xi1+…+βDxiD=x iβ,i=1,2,…,N(3.1)其中x iβ表示向量xi和向量β之间的内积。通常我们会把一个常数1包含在自变量里，以得到简洁表示，也就是说，可以设xi0=1,i=1,2,…,N。对应的β0被称为截距。当D=1，即xi是标量且yi也为简单的标量时，模型称为简单线性回归；当自变量xi是向量时，模型称为多元线性回归。
另外，值得注意的是，自变量不一定是原始的数据特征，可以是原始特征的非线性函数。只要模型关于参数向量β是线性的，模型就被认为是线性模型；如果模型关于参数是非线性的，则被认为是非线性模型。假设(xi)表示对输入特征的变换函数，也称为基函数，那么线性回归可以表示为更一般化的形式，即yi=(xi) β(3.2)常见的基函数有3种，即
① 多项式基函数。j(x)=xj(3.3)② 高斯基函数。j(x)=exp-(x-μj)22s2(3.4)③ S形（sigmoidal）基函数。j(x)=σx-μjs,σ(a)=11+exp(-a)(3.5)3种基函数的示意图如图31所示。
图313种基函数示意图
注: 每张子图展示了使用不同参数产生的3条基函数曲线
得到回归模型后，可以通过公式(3.2)对新的测试数据进行预测,那么测试输入x的预测输出y表示为y=(x) β。为得到该回归函数中参数的最优值，本节将介绍两种对线性回归模型的训练方法: 最小二乘和正则化最小二乘。
【示例】下面通过估计熊猫食量的例子介绍如何使用线性回归建模数据。一篇关于圈养大熊猫食竹量观察的文献记录了4只大熊猫的夜间食竹量，如表31所示。表31大熊猫平均夜间食竹量［1］(单位: kg)熊猫名性别年龄/岁体重1月2月3月4月5月6月7月8月9月10月11月12月莉莉雌10～11102.52.83.32.63.52.74.91.31.71.91.62.53.9青青雌3～482.53.43.73.73.94.15.71.62.12.42.73.34.1金金雄22～23128.01.92.51.72.12.24.51.11.51.21.71.72.1平平雄9～1082.04.24.44.14.64.56.93.23.53.43.43.74.5从表31中可以看出，食竹量与熊猫的性别、年龄、体重及月份都有关系，并且从数据可以简单地分析得到食竹量与月份有两段不同的规律。因此可以以7月为界限，分两段进行线性回归,每一段有24个训练样本。在该示例中，自变量是一个四维向量x=(x1,x2,x3,x4)，因变量是一个标量y。假定用线性模型建模，即yi=β0+β1xi1+β2xi2+β3xi3+β4xi4(3.6)如果要进一步增强模型的灵活性，可以对某些自变量先进行非线性变换，得到新的自变量x~,再进行线性回归建模。此时，虽然模型关于某些变量是非线性的，但是关于参数还是线性的。读者可以尝试对上述示例中的输入特征设置合适的非线性变换，然后进行线性回归建模。
正如前文中介绍的，使用概率模型是建模不确定性的有效方法。概率线性回归的一种实现方式是使用高斯随机噪声实现概率建模。具体来说,观测输出被假设为确定性的线性回归再加上一个高斯随机噪声，表示为y=f(x,β)+，～ (0,σ2)(3.7)其中f(x,β)=(x) β(3.8)根据概率分布的变换关系,可以得到每个观测数据的似然概率分布为p(y|x,β,σ2)= (y|f(x,β),σ2)(3.9)处理实际问题时，模型通常假设数据是独立同分布的，所有观测y的似然概率分布表示为p(y| X,β,σ2)=∏Ni=1 (yi|f(xi,β),σ2)(3.10)确定了模型的概率表示之后，对于新的测试数据，可以使用输出变量的期望作为预测值,计算表达公式为 ［y|x］=∫yp(y|x,β,σ2)dy=f(x,β)(3.11)关于如何确定模型中的参数值，下面将介绍对概率线性回归模型的最大似然估计和最大后验估计，并说明二者分别与最小二乘和正则化最小二乘之间的关系。
3.1.1最小二乘与最大似然
最小二乘法(least square method)中“最小二乘”的意思是最小化误差的平方和，误差是指观测数据的真实输出值和由模型拟合的因变量值之间的差。下面分别给出最小二乘问题的描述，如何求解最小二乘问题，以及概率线性回归的最大似然估计。
(1) 最小二乘问题描述。
给定有N个数据点(xi,yi)的数据集，其中xi为自变量，yi为因变量。模型函数具有形式f(xi,β)，其中β保存了D个可调整的参数。最小二乘问题的目标为调整模型函数的参数最好地拟合数据集。模型对数据的拟合程度是通过其误差来测量的。误差定义为因变量的真实值和模型预测值之间的差，即ei=yi-f(xi,β)(3.12)以曲线拟合为例，误差的几何意义如图32所示。最小二乘法通过最小化平方误差和S学习最优参数值，即S=∑Ni=1e2i=∑Ni=1(yi-f(xi,β))2(3.13)图32误差的几何意义示意图
注: 图中纵向线段长度代表不同数据点的误差
(2) 求解最小二乘问题。
上述平方和的最小化可通过将对优化目标关于参数的导数设为0求解得到。如果分别考虑每一个参数，那么由于模型有D个参数，就有D个梯度方程，即Sβd=0,d=1,2,…,D(3.14)代入公式(3.13)可得-2∑Ni=1yi-f(xi,β)f(xi,β)βd=0,d=1,2,…,D(3.15)公式(3.15)中的梯度方程适用于所有最小二乘问题。每个具体问题有特定的模型表达式和相应的偏导数。当然，求解公式(3.13)给出的最小平方误差和也可以直接使用向量微积分的方法，直接对优化目标关于参数向量求导解得。
下面以线性回归问题为例，具体介绍最小二乘法的解。由公式(3.2)可知，一般化的线性回归模型表示为f(xi,β)=(xi) β。定义X=［x1,x2,…,xN］ ，y=［y1,y2,…,yN］ ，Φ=［(x1),(x2),…,(xN)］ ，那么模型在训练数据上的预测平方误差为S=(y-Φβ) y-Φβ(3.16)根据公式(3.14)可以得到β的最优值满足dSdβ=d((y-Φβ) (y-Φβ))dβ=0 (3.17)其中，0表示元素为0的列向量，d((y-Φβ) (y-Φβ))可利用向量微积分的运算法则（附录C）作进一步化简，即d((y-Φβ) (y-Φβ))=(d(y-Φβ) )(y-Φβ)+(y-Φβ) d(y-Φβ)
=2(y-Φβ) d(y-Φβ)
=-2(y-Φβ) Φdβ
=2(β Φ Φ-y Φ)dβ(3.18)因此，得到β的最优值为β^ s=(Φ Φ)-1Φ y(3.19)确定性线性回归的优化准则通常使用损失函数来定义，最小二乘方法使用的是最小化平方误差和。概率线性回归的优化准则通常以最大似然为目标，这里给出其最大似然解，并说明与最小二乘之间的关系。
(3) 概率线性回归的最大似然估计。
当概率线性回归的似然假设为高斯分布时，如公式(3.10)，其对数似然的表达式可以进一步推导得出lnp(y|X,β,σ2)=-12σ2∑Ni=1yi-f(xi,β)2-N2lnσ2-N2ln(2π)(3.20)最大化公式(3.20)可以获得参数β和σ2的最大似然估计。二者的结果为β^m =(Φ Φ)-1Φ y(3.21)
σ^2m =1N∑Ni=1yi-f(xi,β^m )2(3.22)因此，可以看出当观测数据服从高斯分布时，线性回归参数β的最小二乘解和最大似然估计是等价的。
3.1.2正则化最小二乘与最大后验
回归模型经常遇到数据过拟合（overfitting）问题，也就是模型在训练集上的拟合误差很小，但是在测试集上的误差很大。过拟合通常发生在数据量较少或模型的复杂度太高时。例如，在线性回归中，对数据特征引入多项式变换，回归系数的数量越多，则会导致曲线的波动越大，此时曲线容易对数据产生过拟合。图33给出了4种多项式拟合的效果。
图334种不同的多项式拟合效果
注: 图中小圆圈表示样本，虚线表示真实情况，实线表示拟合曲线，使用的多项式形式为f(x)=∑degj=0wjxj，deg表示多项式的阶数，4张子图分别使用不同的阶数
通常情况下，模型的复杂度是相对的，它与数据量相关。如果训练数据足够多，高复杂度的模型可以很好地拟合数据；如果数据量较少，就需要一个相对简单的模型来拟合数据。在实际应用中，数据量的多少通常是无法改变的，建模者可以控制的是模型的设置，希望可以通过某种约束实现对数据较为合适的拟合。对于使用多项式变换的线性回归的例子，回归系数是关键参数，如果固定多项式的次数，控制回归系数的大小同样可以控制模型复杂度。
(1) 正则化最小二乘。
由于回归系数越大模型波动越大，为了降低过拟合的风险，可以对回归系数进行约束。对最小二乘进行正则化的方法叫作正则化最小二乘。例如，约束回归系数构成的向量的L2范数的平方（‖β‖L2=β β）不超过一个给定值。该约束相当于求解一个带有惩罚项（penalty term）λ‖β‖2的最小二乘的无约束最小化问题。此时，正则化最小二乘的优化目标为S′=∑Ni=1yi-f(xi,β)2+λβ β(3.23)其中λ是常数，可以通过模型选择的方法确定取值。使用L2范数作为惩罚项的正则化最小二乘也称为岭回归［2］。
正则化最小二乘不仅限于L2范数，其他如L1范数（‖β‖L1=∑d|βd|）等也是可行的，不同的正则化项具有不同的约束性质。例如，L2范数的惩罚项可以帮助模型避免过拟合，L1范数的惩罚项除了使得模型减轻过拟合以外，还能够得到较为稀疏的参数解。为了直观理解两种正则化方法，图34展示了这两种正则化项的等高线。
图34L1范数和L2范数的等高线示意图
注: 图中的曲线表示二维空间中的向量x=［x1,x2］ 的L1范数‖x‖L1和L2范数‖x‖L2的等高线
(2) 求解正则化最小二乘问题。
求解正则化最小二乘与求解最小二乘是类似的，同样可以使用求导数的方法得到参数的闭式解。对于使用L2范数的正则化最小二乘，其最优解满足dS′dβ=dy-Φβ y-Φβ+λβ βdβ=0 (3.24)d（（y-Φβ） （y-Φβ））可利用向量微积分的运算法则（附录C）进一步化简为d((y-Φβ) (y-Φβ)+λβ β)=2((y-Φβ) d(y-Φβ)+λβ dβ)
=-2((y-Φβ) Φ-λβ )dβ
=2(β Φ Φ-y Φ+λβ )dβ(3.25)因此，得到β的最优值为β^rls=(λI+Φ Φ)-1Φ y(3.26)(3) 概率线性回归的最大后验估计。
回顾3.1.1节，线性回归的似然假设为高斯分布时，如果使用最大后验估计来获取模型参数，需要假设参数的先验分布。在高斯似然的模型中，通常使用高斯分布作为先验，这样得到的概率线性回归中参数的后验分布还是高斯分布。一种简单常用的先验分布为p(β|α)= (β|0,α-1I)(3.27)根据贝叶斯公式可以得出参数的对数后验分布为lnp(β|X,y,α,σ2)=-12σ2∑Ni=1yi-f(xi,β)2-α2β β+const(3.28)其中const表示与β无关的项。将最大后验估计的优化目标(公式(3.28))与正则化最小二乘的目标(公式(3.23))对比可以发现，当α=λ/σ2时，二者的解是等价的。因此，通过给模型参数增加先验并进行最大后验估计的方法同样可以达到减轻过拟合的效果。
3.2贝叶斯线性回归
3.1节介绍的使用最大似然估计与最大后验估计的线性回归模型，模型参数的值完全通过数据集训练得出。一旦得到β^，可以根据模型估计任意新数据点x的输出值: y^=(x) β^。它们得到的是参数的点估计，即给定数据时可能性最大的估计。但是，当数据集比较小或不确定性较大时，将估计表示为一个可能值的分布更加合理。这样得到的输出预测值将是一个分布。下面介绍可以给出参数与预测值分布的贝叶斯线性回归。
贝叶斯线性回归（Bayesian linear regression）将贝叶斯框架应用到线性回归中，回归系数β被假设为有一特定先验分布的随机变量，此先验分布可以影响回归系数的解。另外，贝叶斯参数估计不是给出回归系数的最佳单点估计，而是给出完整的后验分布，这种方式描述了估计量的不确定性。
考虑一个标准的线性回归问题，对于i=1,2,…,N, 假设在给定自变量xi的情况下,yi产生的公式为yi=x iβ+i(3.29)其中β是D×1维向量，i是独立同分布的随机变量，并且i～ (0,σ2)。定义X=［x1,x2,…,xN］ ，y=［y1,y2,…,yN］ ，可以得到因变量y的似然函数为
p(y|X,β,σ2)∝(σ2)-N2exp-12σ2(y-Xβ) (y-Xβ)(3.30)
即y～ (Xβ,σ2I)。
在贝叶斯方法中，参数的先验概率分布为模型提供了额外信息。先验可以根据领域知识和已知信息采取不同的函数形式。确定似然函数后，对于一个任意的先验分布，后验分布不一定存在解析形式。这里讨论可以使后验分布被解析地推导出来的情况，常用的方法是设置似然函数的共轭先验。下面给出共轭先验的定义。
如果先验分布和似然函数可以使后验分布和先验分布具有相同的形式，就称先验分布与似然函数是共轭的，该先验称为该似然函数的共轭先验。共轭的好处是让后验分布与先验分布具有相同的形式，从而便于求解。
以上文介绍的线性回归为例，给定模型的似然假设公式(3.30)，需要进行贝叶斯估计的参数包括β和σ2。为了使得后验分布可以得到与先验分布相同的形式，这里假设参数β和σ2的联合先验为p(β,σ2)=p(σ2)p(β|σ2)(3.31)其中p(σ2)是逆伽马分布InvGamma(a0,b0)，即p(σ2)∝(σ2)-a0-1exp-b0σ2(3.32)而p(β|σ2)的条件先验密度服从正态分布 (μ0,σ2Λ-10)，即p(β|σ2)∝exp-12σ2(β-μ0) Λ0(β-μ0)(3.33)给定β和σ2的先验假设，根据贝叶斯公式，可以得到贝叶斯线性回归参数的后验分布为p(β,σ2|y,X)=p(β|σ2,y,X)p(σ2|y,X)∝
p(y|X,β,σ2)p(β|σ2)p(σ2)(3.34)将公式(3.30)、(3.32)和(3.33)代入(3.34)，可得p(β|σ2,y,X)是高斯分布 (β|μN,σ2Λ-1N)，以及p(σ2|y,X)是逆伽马分布InvGamma(σ2|aN,bN)，其参数的具体表示为ΛN=(X X+Λ0)
μN=(ΛN)-1(X y+Λ0μ0)
aN=a0+N2
bN=b0+12(y y+μ 0Λ0μ0-μ NΛNμN)(3.35)3.3逻 辑 回 归
前文介绍了使用线性函数预测连续取值的变量，这类问题称为回归问题。很多时候，也需要预测离散取值变量，例如判断一张图像属于哪个目标类别，这变成了分类问题。逻辑回归在线性回归的基础上实现了二类和多类分类。
逻辑回归（logistic regression）［3］模型是一种常用的分类算法，也可以认为是一种因变量为离散值的回归模型。逻辑回归可以处理二类分类和多类分类问题。在二类逻辑回归中，因变量只有两种取值，例如“0”或“1”。在多类逻辑回归中，因变量有两种以上的离散取值。
本节首先介绍二类逻辑回归，然后介绍多类逻辑回归。
3.3.1二类逻辑回归
二类逻辑回归模型使用一个或多个自变量（特征）来估计因变量取值的概率。输出通常被编码为“0”或“1”。模型本身根据输入仅仅建模了输出的概率，并不执行分类，即模型本身并不是一个分类器。当然，通常可以使用此模型构造一个分类器，例如，选择一个阈值，将概率大于此阈值的输入分为一类，小于此阈值的分为另一类。逻辑回归模型使用逻辑函数（logistic function），将线性回归的返回值转换为区间［0,1］内的值，用于表示自变量属于某个类别的概率，即因变量取值为“0”或“1”的概率。
逻辑函数也称为sigmoid函数，输入可以是任意实数x(x∈R)，输出的值属于区间［0,1］。逻辑函数σ(x)的表达式为σ(x)=exex+1=11+e-x(3.36)其函数曲线如图35所示。它是一个S形曲线，在横坐标取值远离0时，纵坐标的值趋近0或1。
图35逻辑函数示意图
逻辑回归使用逻辑函数和回归模型可以解决二分类问题，其中逻辑函数的返回值用于表示二分类问题中的正类或负类的概率。假设f是自变量x的一个线性函数，即f=θ x。逻辑回归假设样本x属于正类的概率为p(y=1|x)=hθ(x)=σ(θ x)=11+exp(-θ x)(3.37)那么，x属于负类的概率为
p(y=0|x)=1-p(y=1|x)=1-hθ(x)=11+exp(θ x)(3.38)
逻辑回归可以从两个角度定义目标函数，一种是从最大似然的角度，一种是从直接构建损失的角度。逻辑回归的优化目标是学习得到合适的参数值，使得概率p(y=1|x)=hθ(x)在当x属于“1”类时值比较大，且p(y=0|x)=1-hθ(x)在当x属于“0”类时值比较大。
从最大似然的角度分析，假设每一个样本的类标签都是独立同分布的伯努利变量，伯努利变量取值为1和0的概率分别为公式(3.37)和公式(3.38)。对于有标签的训练集{(xi,yi): i=1,2,…,N}，N个独立样本的联合似然可以写成p(y|θ)=∏Ni=1p(yi=1|xi)yi(1-p(yi=1|xi))(1-yi)(3.39)最大化似然等价于最小化负对数似然，因此，最大似然得到的损失函数为-lnp(y|θ)=-∑Ni=1［yilnp(yi=1|xi)+
(1-yi)ln(1-p(yi=1|xi))］(3.40)从构建损失函数的角度分析，逻辑回归使用真实概率分布与模型概率分布的交叉熵损失来直接定义训练集{(xi,yi): i=1,2,…,N}的损失函数。假设每个样本的真实分布为q(yi|xi),那么，q(yi=1|xi)=yi，且q(yi=0|xi)=1-yi。分布q(yi|xi)和p(yi|xi)的交叉熵为H(q(yi|xi),p(yi|xi))=-∑yiq(yi|xi)lnp(yi|xi)(3.41)因此,逻辑回归的交叉熵损失为J(θ)=∑Ni=1H［q(yi|xi),p(yi|xi)］
=-∑Ni=1yilnhθ(xi)+(1-yi)ln(1-hθ(xi))(3.42)无论从最大似然角度还是最小损失函数角度，二者得到的目标损失是一致的。可以通过最小化J(θ)找到假设函数hθ(x)中θ的最优值，从而学得分类器。关于该目标的优化得不到闭式解［4］，因此常用基于梯度的迭代优化方法，例如一阶梯度下降或基于二阶梯度的牛顿法等。使用梯度下降等方法优化θ，需要计算J(θ)关于θ的梯度，计算公式为  θJ(θ)=dJ(θ)dθ =∑Ni=1［(σ(θ xi)-yi)x i］dθdθ 
=∑ixi［hθ(xi)-yi］(3.43)其中，梯度的运算过程利用了sigmoid函数的导数性质: dσ(x)=σ(x)(1-σ(x))dx。使用牛顿法进行优化,则还需要计算Hessian矩阵。
得到合适的参数值后，对于新的测试样本x，如果p(y=1|x)>p(y=0|x)，那么将此样本标记为“1”类，否则标记为“0”类。相应的决策函数为: 如果p(y=1|x)>0.5，那么y=1。通常情况下,选择0.5作为阈值进行决策，在很多实际应用中，也可以根据特定的情况选择不同的阈值。例如，如果对正例的判别查准率要求高，可以选择大于0.5的值作为阈值；如果对正例的查全率要求高，可以选择小于0.5的值作为阈值。
3.3.2多类逻辑回归
多类逻辑回归（multinomial logistic regression）的基本原理与二类逻辑回归类似，差别在于多类逻辑回归中因变量yi的取值可以大于两个，一个C类逻辑回归的因变量可以在1～C取任意一个整数。多类逻辑回归使用softmax实现从实数到类别概率的转换。
定义类别标签为c∈{1,2,…,C}，每一个类别对应于一个回归函数,即fc(xi)=θ cxi(3.44)其中θc是与类别c对应的回归系数，xi是第i个样本向量。经过softmax函数转换后得到样本属于某一类别的概率为p(yi=c)=exp(θ cxi)∑Ck=1exp(θ kxi)(3.45)根据公式（3.45），样本被分为概率最大的那一类。每个向量θc中未知的参数可以通过最大似然或最小化交叉熵进行优化。多类逻辑回归的似然函数为p(y|θ1,θ2,…,θC)=∏Ni=1∏Cc=1p(yi=c|xi)I(yi=c)(3.46)其中，I(yi=c)仅当yi=c时函数值为1，其余为0。对应的负对数似然，也就是交叉熵损失为
-lnp(y|θ1,θ2,…,θK)=-∑Ni=1∑Cc=1I(yi=c)lnp(yi=c|xi)(3.47)
与二类逻辑回归类似，由于优化目标中包含非线性函数，通常得不到闭式解，因此常用的方法是基于梯度的迭代优化。此外，无论是二类逻辑回归还是多类逻辑回归，使用最大后验估计或最小化带惩罚项的交叉熵损失可以防止模型过拟合。
3.4贝叶斯逻辑回归
本节以两分类为例介绍贝叶斯逻辑回归（Bayesian logistic regression）。逻辑回归是一种判别式概率线性分类器p(y=1|x,θ)=σ(θ x)。贝叶斯逻辑回归通过贝叶斯参数估计学习参数的后验分布，并且利用该分布进行预测。
已知观测数据X=［x1,x2,…,xN］ ，y=［y1,y2,…,yN］ ，逻辑回归使用的似然分布导致后验分布p(θ|X,y)难以有解析表达，因此通常使用其他典型分布q(θ)来近似后验分布。预测时，即便使用了近似分布，对新样本x的预测分布p(y=1|x)≈∫σ(θ x)q(θ)dθ的估计仍然是难解的。因此，贝叶斯逻辑回归通常使用近似求解方法。
一方面，后验分布p(θ|X,y)等于先验乘以似然,再进行归一化。其中先验通常假设为p(θ)= (θ|m0,S0)(3.48)逻辑回归的似然为p(y|X,θ)=∏Ni=1p(yi=1|xi)yi(1-p(yi=1|xi))1-yi(3.49)对逻辑回归中的后验分布进行精确求解非常困难，这时可以通过使用拉普拉斯近似得到近似的高斯后验分布q(θ)。
另一方面，预测分布p(y=1|x)≈∫σ(θ x)q(θ)dθ 需要关于sigmoid函数和高斯分布的乘积求积分，其精确求解也是十分困难的，可通过将sigmoid函数用逆probit函数近似得到其近似解［4］。下面对这两方面的近似进行详细介绍。
（1） 拉普拉斯近似。
对后验分布的拉普拉斯近似是通过数值优化算法得到一个以θ0为均值的高斯分布q(θ)，作为真实后验的近似分布为q(θ)=1(2π)D/2|SN|1/2exp-12(θ-θ0) S-1N(θ-θ0)= (θ|θ0,SN)(3.50)其中，均值θ0是真实后验分布的最大值对应的变量值，协方差矩阵是负对数真实后验分布-lnp(θ|X,y)的Hessian矩阵（附录C）在θ=θ0处的逆，即SN=-  lnp(θ|X,y)|θ=θ0-1。下面来看均值θ0和协方差矩阵SN的具体计算过程。
已知参数服从高斯先验p(θ)= (θ|m0,S0)，其中m0和S0是超参数。后验分布p(θ|X,y)∝p(θ)p(y|X,θ)。将先验概率（公式(3.48)）和逻辑回归的似然函数（公式(3.49)）代入贝叶斯公式可得lnp(θ|X,y)=-12(θ-m0) S-10(θ-m0)+∑Ni=1［yilnp(yi=1|xi,θ)+
(1-yi)ln(1-p(yi=1|xi,θ))］+const(3.51)最大化该对数后验分布lnp(θ|X,y),可以得到参数的最大后验估计θmap，作为近似分布q(θ)的均值。-lnp(θ|X,y)的Hessian矩阵计算为H=-  lnp(θ|X,y)=d2lnp(θ|X,y)dθdθ 
=dTr［(θ-m0) S-10dθ］-d∑Ni=1(yi-σ(θ xi))x idθdθdθ 
=Tr［S-10dθdθ ］+Tr∑Ni=1σ(θ xi)(1-σ(θ xi))xix idθdθ dθdθ 
=S-10+∑Ni=1p(yi=1|xi,θ)(1-p(yi=1|xi,θ))xix i(3.52)其中运算过程利用了sigmoid函数的导数性质: dσ(x)=σ(x)1-σ(x)dx。得到H之后，根据SN=(H|θ=θmap)-1得到近似分布的协方差矩阵SN，可以得到后验分布的高斯近似q(θ)= (θ|θmap,SN)。
(2) 逆probit函数近似。
得到近似后验分布后，对于给定的新特征向量x，其属于类别“1”的预测分布可以通过似然关于后验p(θ|X,y)的积分得到，即p(y=1|x)=∫p(y=1,θ|x)dθ
=∫p(y=1|x,θ)p(θ|X,y)dθ
≈∫σ(θ x)q(θ)dθ(3.53)属于类别“0”的概率为p(y=0|x)=1-p(y=1|x)(3.54)下面对公式(3.53)作进一步化简，由于函数σ(θ x)仅通过θ x的值依赖于θ，因此定义新的变量a=θ x，并引入Dirac delta函数δ(·)，可以得到σ(θ x)≈∫δ(a-θ x)σ(a)da。因此，公式(3.53)的结果可以表示为∫σ(θ x)q(θ)dθ=∫∫δ(a-θ x)σ(a)daq(θ)dθ
=∫σ(a)∫δ(a-θ x)q(θ)dθda(3.55)其中，∫δ(a-θ x)q(θ)dθ是关于a的函数，并且可验证为是一个高斯概率分布，记为p(a)= (a|μa,σ2a)，其中均值与方差分别为μa= ［a］=∫p(a)ada=∫q(θ)θ xdθ=θ mapx(3.56)
σ2a=var［a］=∫p(a)a2da- ［a］2=∫q(θ)(θ x)2dθ-(θ mapx)2
=x SNx(3.57)预测分布可以表示为p(y=1|x)=∫σ(a)p(a)da=∫σ(a) (a|μa,σ2a)da(3.58)注意,在公式(3.58)的积分中，关于sigmoid和Gaussian乘积的积分是不可解的，通常使用逆probit函数来替代sigmoid函数。定义标准高斯分布的累积分布函数为Φ(a)=∫a-∞ (w|0,1)dw(3.59)该函数也称为逆probit函数。由于累积分布函数的值域是（0,1），因此可用逆probit函数来近似sigmoid函数。为使二者尽可能一致，需对逆probit函数的自变量进行放缩，即使用Φ(λa)来近似σ(a)，并且λ通常设置为λ=π/8，此时两者在原点具有相同的斜率（即导数相同）［4］。高斯分布和逆probit函数相乘后的积分还是一个逆probit函数，即∫Φ(λa) (a|μa,σ2a)da=Φμa(λ-2+σ2a)1/2(3.60)将公式（3.60)应用到公式(3.58)中,可以获得最终的预测概率为