MLAPP 读书笔记 - 08 逻辑回归(Logistic regression)
A Chinese Notes of MLAPP,MLAPP 中文笔记项目
https://zhuanlan.zhihu.com/python-kivy
转载：cycleuser
侵权立删，本意是为广大学习爱好者提供中文资料，同时在这里补充了图片，方便阅读，如果有任何法律问题，我会立即删除。
8.1 概论
构建概率分类器有一种方法是建立形式为 p ( y , x ) p(y,x) p(y,x)的联合模型,然后以x为条件,推 p ( y ∣ x ) p(y|x) p(y∣x).这叫生成方法(generative approach).另外一种办法是直接以 p ( y ∣ x ) p(y|x) p(y∣x)的形式去拟合一个模型.这就叫做启发式方法(discriminative approach),本章就要讲这个方法.具体来说就是要假设有一些参数为线性的启发模型.我们会发现这样的模型拟合起来特别简单.在8.6,我们会对生成方法和启发方法进行对比,在本书后面的章节中,还会讲到非线性的启发模型和非参数化的启发模型.
8.2 模型选择
正如在本书1.4.6讲过的,逻辑回归对应的是下面这种二值化分类模型:
p
(
y
∣
x
,
w
)
=
B
e
r
(
y
∣
s
i
g
m
(
w
T
x
)
)
p(y|x,w)=Ber(y|sigm(w^Tx))
p(y∣x,w)=Ber(y∣sigm(wTx))(8.1)
图1.19(b)所示的就是一个一维例子.逻辑回归可以很容易扩展用于高维度输入特征向量.比如图8.1所示的就是对二维输入和不同的权重向量w的逻辑回归 p ( y = 1 ∣ x , w ) = s i g m ( w T x ) p(y=1|x,w)=sigm(w^Tx) p(y=1∣x,w)=sigm(wTx).如果设置概率0.5的位置为阈值,就能得到一个决策边界,其范数(垂线)为w.
8.3 模型拟合
本节讲对逻辑回归模型进行参数估计的算法.
此处参考原书图8.1
8.3.1 最大似然估计(MLE)
逻辑回归的负对数似然函数为:
N
L
L
(
w
)
=
−
∑
i
=
1
N
log
[
μ
i
I
(
y
i
=
1
)
×
(
1
−
μ
i
)
I
(
y
i
=
0
)
]
(8.2)
=
−
∑
i
=
1
N
log
[
y
i
log
μ
i
+
(
1
−
y
i
)
log
(
1
−
μ
i
)
]
(8.3)
\begin{aligned} NLL(w)&= -\sum^N_{i=1}\log [\mu_i^{I(y_i=1)}\times (1-\mu_i)^{I(y_i=0)}] &\text{(8.2)}\\ &= -\sum^N_{i=1}\log [y_i\log \mu_i+(1-y_i)\log(1-\mu_i)] &\text{(8.3)}\\ \end{aligned}
NLL(w)=−i=1∑Nlog[μiI(yi=1)×(1−μi)I(yi=0)]=−i=1∑Nlog[yilogμi+(1−yi)log(1−μi)](8.2)(8.3)
这也叫交叉熵误差函数(cross-entropy error function)参考本书2.8.2.
这个公式还有另外一个写法,如下所示.设 y ~ i ∈ { − 1 , + 1 } \tilde y_i \in \{-1,+1\} y~i∈{−1,+1},而不是 y i ∈ { 0 , 1 } y_i\in\{0,1\} yi∈{0,1}.另外还有 p ( y = − 1 ) = 1 1 + exp ( − w T x ) p(y=-1)=\frac{1}{ 1+\exp (-w^Tx)} p(y=−1)=1+exp(−wTx)1和 p ( y = 1 ) = 1 1 + exp ( + w T x ) p(y=1)=\frac{1}{ 1+\exp (+w^Tx)} p(y=1)=1+exp(+wTx)1.这样有:
N L L ( w ) = ∑ i = 1 N log ( 1 + exp ( − y ~ i w T x i ) ) NLL(w)=\sum^N_{i=1}\log(1+\exp(-\tilde y_i w^Tx_i)) NLL(w)=∑i=1Nlog(1+exp(−y~iwTxi))(8.4)
和线性回归不一样的是,在逻辑回归里面,我们不再能以闭合形式写出最大似然估计(MLE).所以需要使用一个优化算法来计算出来.为了这个目的,就要推到梯度(gradient)和海森矩阵(Hessian).
此处参考原书图8.2
如练习8.3所示,很明显梯度和海森矩阵分别如下所示:
g
=
d
d
w
f
(
w
)
=
∑
i
∗
μ
i
−
y
i
)
x
i
=
X
T
(
μ
−
y
)
(8.5)
H
=
d
d
w
g
(
w
)
T
=
∑
i
(
∇
w
μ
i
)
x
i
T
=
∑
i
μ
i
(
1
−
μ
i
)
x
x
x
i
T
(8.6)
=
X
T
S
X
(8.7)
\begin{aligned} g &=\frac{d}{dw}f(w)=\sum_i*\mu_i-y_i)x_i=X^T(\mu-y) &\text{(8.5)}\\ H &= \frac{d}{dw}g(w)^T=\sum_i(\nabla_w\mu_i)x_i^T=\sum_i\mu_i(1-\mu_i)x_xx_i^T &\text{(8.6)}\\ &= X^TSX &\text{(8.7)}\\ \end{aligned}
gH=dwdf(w)=i∑∗μi−yi)xi=XT(μ−y)=dwdg(w)T=i∑(∇wμi)xiT=i∑μi(1−μi)xxxiT=XTSX(8.5)(8.6)(8.7)
其中的 S = △ d i a g ( μ i ( 1 − μ i ) ) S\overset{\triangle}{=}diag(\mu_i(1-\mu_i)) S=△diag(μi(1−μi)).通过练习8.3也可以证明海森矩阵H是正定的(positive definite).因此负对数似然函数NLL就是凸函数(convex)有唯一的全局最小值.接下来就说一下找到这个最小值的方法.
8.3.2 梯度下降(gradient descent)
无约束优化问题的最简单算法,可能就是梯度下降(gradient descent)了,也叫做最陡下降(Steepest descent).写作下面的形式:
θ k + 1 = θ k − η k g k \theta_{k+1}=\theta_k -\eta_kg_k θk+1=θk−ηkgk(8.8)
其中的 η k \eta_k ηk是步长规模(step size)或者也叫学习率(learning rate).梯度下降法的主要问题就是:如何设置步长.这个问题还挺麻烦的.如果使用一个固定的学习率,但又太小了,那收敛就会很慢,但如果要弄太大了呢,又可能最终不能收敛.这个过程如图8.2所示,其中对下面的(凸函数)进行了投图:
f ( θ ) = 0.5 ( θ 1 2 − θ 2 ) + 0.5 ( θ 1 ) 2 f(\theta)=0.5(\theta_1^2-\theta_2)+0.5(\theta_1)^2 f(θ)=0.5(θ12−θ2)+0.5(θ1)2(8.9)
任意设置从(0,0)开始.在图8.2(a)中,使用了固定步长为 η = 0.1 \eta =0.1 η=0.1;可见沿着低谷部位移动很慢.在图8.2(b)中,使用的是固定步长 η = 0.6 \eta =0.6 η=0.6;很明显这样一来算法很快就跑偏了,根本就不能收敛了.
此处参考原书图8.3
所以就得想个靠谱的办法来选择步长,这样才能保证无论起点在哪里最终都能收敛到局部最优值.(这个性质叫做全局收敛性(global convergence),可千万别跟收敛到全局最优值弄混淆哈.)通过泰勒定理(Taylor’s theorem),就得到了:
f ( θ + η d ) ≈ f ( θ ) + η g T d f(\theta+\eta d)\approx f(\theta)+\eta g^Td f(θ+ηd)≈f(θ)+ηgTd(8.10)
其中的d是下降方向.所以如果 η \eta η足够小了,则 f ( θ + η d ) < f ( θ ) f(\theta+\eta d)< f(\theta) f(θ+ηd)<f(θ),因为梯度会是负值的.不过我们并不希望步长太小,否则就要运算很久才能到达最小值了.所以要选一个能够最小化下面这个项的步长 η \eta η:
ϕ ( η ) = f ( θ k + η d k ) \phi(\eta)=f(\theta_k +\eta d_k) ϕ(η)=f(θk+ηdk)(8.11)
这就叫线性最小化(line minimization)或者线性搜索(line search).有很多种方法来借这个一维优化问题,具体细节可以参考(Nocedal and Wright 2006).
图8.3(a)展示了上面那个简单问题中的线性搜索.不过我们会发现线性搜索得到的梯度下降路径会有一种扭折行为(zig-zag behavior).可以看到其中一次特定线性搜索要满足 η k = arg min η > 0 ϕ ( η ) \eta_k =\arg \min_{\eta>0}\phi (\eta) ηk=argminη>0ϕ(η).优化的一个必要条件就是导数为零,即 ϕ ′ ( η ) = 0 \phi'(\eta)=0 ϕ′(η)=0.通过链式规则(chain rule), ϕ ′ ( η ) = d T g \phi'(\eta)=d^Tg ϕ′(η)=dTg,其中 g = f ′ ( θ + η d ) g=f'(\theta+\eta d) g=f′(θ+ηd)是最后一步的梯度.所以要么就有 g = 0 g=0 g=0,意思就是已经找到了一个固定点(stationary point);要么就是$g\perp d $,意味着这一步所在点的位置上局部梯度和搜索方向相互垂直.因此连续起来方向就是正交的,如图8.3(b)所示,这就解释了搜索路径的扭折行为.
降低这种扭折效应的一种简单的启发式方法就是增加一个动量项(momentum term) θ k − θ k − 1 \theta_k-\theta_{k-1} θk−θk−1:
θ k + 1 = θ k − η k g k + μ k ( θ k − θ k − 1 ) \theta_{k+1}=\theta_k-\eta_kg_k+\mu_k(\theta_k-\theta_{k-1}) θk+1=θk−ηkgk+μk(θk−θk−1)(8.12)
上式中的 0 ≤ μ k ≤ 1 0\le \mu_k \le 1 0≤μk≤1控制了动量项的重要程度.在优化领域中,这个方法叫做重球法(heavy ball method ,参考 Bertsekas 1999).
另外一种最小化扭折行为的方法是使用共轭梯度(conjugate gradients)(参考Nocedal and Wright 2006第五章,或者Golub and van Loan 1996,10.2).这个方法是选择形式为 f ( θ ) = θ T A θ f(\theta)=\theta^TA\theta f(θ)=θTAθ的二次形式(quadratic objectives),这是在解线性方程组的时候出现的.不过非线性的共轭梯度就不太受欢迎了.
8.3.3 牛顿法
算法8.1 最小化一个严格凸函数的牛顿法
- 初始化一个 θ 0 \theta_0 θ0;
- 对于k=1,2,…等,一直到收敛为止,重复下面步骤:
- 估算 g k = ∇ f ( θ k ) g_k=\nabla f(\theta_k) gk=∇f(θk)
- 估算 H k = ∇ 2 f ( θ k ) H_k=\nabla^2 f(\theta_k) Hk=∇2f(θk)
- 对 d k d_k dk求解 H k d k = − g k H_kd_k=-g_k Hkdk=−gk
- 使用线性搜索来找到沿着 d k d_k dk方向的步长 η k \eta_k ηk
- θ K + 1 = θ k + η k d k \theta_{K+1}=\theta_k+\eta_kd_k θK+1=θk+ηkdk
如果把空间曲率(curvature)比如海森矩阵(Hessian)考虑进去,可以推导出更快速的优化方法.这样方法就成了二阶优化方法了(second order optimization metods).如果不考虑曲率的那个就叫做牛顿法(Newton’s algorithm).这是一个迭代算法(iterative algorithm),其中包含了下面形式的更新步骤:
θ
k
+
1
=
θ
k
−
η
k
H
k
−
1
g
k
\theta_{k+1}=\theta_k-\eta_kH_k^{-1}g_k
θk+1=θk−ηkHk−1gk(8.13)
完整的伪代码如本书算法2所示.
这个算法可以按照下面步骤推导.设构建一个二阶泰勒展开序列来在 θ k \theta_k θk附近估计 f ( θ ) f(\theta) f(θ):
f q u a d ( θ ) = f k + g k T ( θ − θ k ) + 1 2 ( θ − θ k ) T H k ( θ − θ k ) f_{quad}(\theta)=f_k+g^T_k(\theta-\theta_k)+\frac{1}{2}(\theta-\theta_k)^TH_k(\theta-\theta_k) fquad(θ)=fk+gkT(θ−θk)+21(θ−θk)THk(θ−θk)(8.14)
重写成下面的形式
f q u a d ( θ ) = θ T A θ + b T θ + c f_{quad}(\theta)=\theta^TA\theta+b^T\theta+c fquad(θ)=θTAθ+bTθ+c(8.15)
其中:
A = 1 2 H k , b = g k − H k θ k , c = f k − g k T θ k + 1 23 θ k T H k θ k A=\frac{1}{2}H_k,b=g_k-H_k\theta_k,c=f_k-g^T_k\theta_k+\frac{1}{23}\theta^T_kH_k\theta_k A=21Hk,b=gk−Hkθk,c=fk−gkTθk+231θkTHkθk(8.16)
f q u a d f_{quad} fquad最小值为止在:
θ = − 1 2 A − 1 b = θ k − H k − 1 g k \theta=-\frac{1}{2}A^{-1}b=\theta_k-H_k^{-1}g_k θ=−21A−1b=θk−Hk−1gk(8.17)
因此牛顿步长 d k = − H k − 1 g k d_k=-H_k^{-1}g_k dk=−Hk−1gk就可以用来加到 θ k \theta_k θk上来最小化在 θ k \theta_k θk附近对 f f f的二阶近似.如图8.4(a)所示.
此处参考原书图8.4
在最简单的形式下,牛顿法需要海森矩阵 H k H_k Hk为正定矩阵,这保证了函数是严格凸函数.否则,目标函数非凸函数,那么海森矩阵 H k H_k Hk就可能不正定了,所以 d k = − H k − 1 g k d_k=-H_k^{-1}g_k dk=−Hk−1gk就可能不是一个下降方向了(如图8.4(b)所示).这种情况下,简单的办法就是逆转最陡下降方向, d k = − g k d_k=-g_k dk=−gk.列文伯格-马夸特算法(Levenberg Marquardt algorithm)是一种在牛顿步长和最陡下降步长之间这种的自适应方法.这种方法广泛用于解非线性最小二乘问题.一种替代方法就是:不去直接计算 d k = − H k − 1 g k d_k=-H_k^{-1}g_k dk=−Hk−1gk,可以使用共轭梯度来解关于 d k d_k dk的线性方程组 H k d k = − g k H_kd_k=-g_k Hkdk=−gk.如果 H k H_k Hk不是正定矩阵,只要探测到了负曲率,就可以简单地截断共轭梯度迭代,这样就叫做截断牛顿法(truncated Newton).
8.3.4 迭代重加权最小二乘法(Iteratively reweighted least squares,缩写为IRLS)
接下来试试将牛顿法用到二值化逻辑回归中求最大似然估计(MLE)上面.在这个模型中第 k + 1 k+1 k+1次迭代中牛顿法更新如下所示(设 η k = 1 \eta_k=1 ηk=1,因此海森矩阵(Hessian)是确定的):
w k + 1 = w k − H − 1 g k (8.18) = w k + ( X T S k X ) − 1 X T ( y − μ k ) (8.19) = ( X T S k X ) − 1 [ ( X T S k X ) w k + X T ( y − μ k ) ] (8.20) = ( X T S k X ) − 1 X T [ S k X w k + y − μ k ] (8.21) = ( X T S k X ) − 1 X T S k z k (8.22) \begin{aligned} w_{k+1}&= w_k-H^{-1}g_k &\text{(8.18)}\\ &= w_k+(X^TS_kX)^{-1}X^T(y-\mu_k) &\text{(8.19)}\\ &= (X^TS_kX)^{-1}[(X^TS_kX)w_k+X^T(y-\mu_k)] &\text{(8.20)}\\ &= (X^TS_kX)^{-1}X^T[S_kXw_k+y-\mu_k] &\text{(8.21)}\\ &= (X^TS_kX)^{-1}X^TS_kz_k &\text{(8.22)}\\ \end{aligned} wk+1=wk−H−1gk=wk+(XTSkX)−1XT(y−μk)=(XTSkX)−1[(XTSkX)wk+XT(y−μk)]=(XTSkX)−1XT[SkXwk+y−μk]=(XTSkX)−1XTSkzk(8.18)(8.19)(8.20)(8.21)(8.22)
然后就可以定义工作响应函数(working response)如下所示:
z
k
=
△
X
w
k
+
S
k
−
1
(
y
−
μ
k
)
z_k\overset{\triangle}{=} Xw_k+S_k^{-1}(y-\mu_k)
zk=△Xwk+Sk−1(y−μk)(8.23)
等式8.22就是一个加权最小二乘问题(weighted least squares problem),是要对下面的项最小化:
∑ i = 1 N S k i ( z k i − w T x i ) 2 \sum^N_{i=1}S_{ki}(z_{ki}-w^Tx_i)^2 ∑i=1NSki(zki−wTxi)2(8.24)
由于 S k S_k Sk是一个对角矩阵,所以可以把目标函数写成成分的形式(对每个 i = 1 : N i=1:N i=1:N):
z k i = w k T x i + y i − μ k i μ k i ( 1 − μ k i ) z_{ki}=w_k^Tx_i+\frac{y_i-\mu_{ki}}{\mu_{ki}(1-\mu_{ki})} zki=wkTxi+μki(1−μki)yi−μki(8.25)
这个算法就叫做迭代重加权最小二乘法(Iteratively reweighted least squares,缩写为IRLS),因为每次迭代都解一次加权最小二乘法,其中的权重矩阵 S k S_k Sk在每次迭代都变化.伪代码参考本书配套算法10.
算法8.2 迭代重加权最小二乘法(IRLS)
- w = 0 D w=0_D w=0D;
- w 0 = log ( y ˉ / ( 1 − y ˉ ) ) w_0=\log(\bar y/ (1-\bar y)) w0=log(yˉ/(1−yˉ))
- 重复下面步骤:
-
$\eta_i=w_0+w^Tx_i$ -
$\mu_i=sigm(\eta_i)$ -
$s_i=\mu_i(1-\mu_i)$ -
$z_i=\eta_i+\frac{y_i-\mu_i}{s_i}$ -
$S=diag(s_{1:N})$ -
$w=(X^TSX)^{-1}X^TSz$ - 直到收敛
8.3.5 拟牛顿法
二阶优化算法的源头都是牛顿法,在本书8.3.3中讲到过.不过很不幸的是计算出来海森矩阵H的运算开销成本太高了.拟牛顿法(Quasi-Newton methods)就应运而生了,以迭代方式使用从每一步的梯度向量中学到的信息来构建对海森矩阵的估计.最常用的方法就是BFGS方法(这四个字母是发明这个算法的四个人的名字的首字母Broyden, Fletcher, Goldfarb, Shanno),这个方法是使用下面所示定义的
B
k
≈
H
k
B_k\approx H_k
Bk≈Hk来对海森矩阵进行估计:
B
k
+
1
=
B
k
+
y
k
y
k
T
y
k
T
s
k
−
(
B
k
s
k
)
(
B
k
s
k
)
T
s
k
T
B
k
s
k
(8.26)
s
k
=
θ
k
−
θ
k
−
1
(8.27)
y
k
=
g
k
−
g
k
−
1
(8.28)
\begin{aligned} B_{k+1}& =B_k+\frac{y_ky_k^T}{y_k^Ts_k}-\frac{(B_ks_k)(B_ks_k)^T}{s_k^TB_ks_k} &\text{(8.26)}\\ s_k& = \theta_k-\theta_{k-1} &\text{(8.27)}\\ y_k& = g_k-g_{k-1} &\text{(8.28)}\\ \end{aligned}
Bk+1skyk=Bk+ykTskykykT−skTBksk(Bksk)(Bksk)T=θk−θk−1=gk−gk−1(8.26)(8.27)(8.28)
这是对矩阵的二阶更新(rank-two update),这确保了矩阵保持正定(在每个步长的特定限制下).通常使用一个对角线估计来启动算法,即设 B 0 = I B_0=I B0=I.所以BFGS方法可以看做是对海森矩阵使用对角线加上低阶估计的方法.
另外BFGS方法也可以对海森矩阵的逆矩阵进行近似,通过迭代更新 C l ≈ = H k − 1 C_l\approx = H_k^{-1} Cl≈=Hk−1,如下所示:
C k + 1 = ( I − s k s k T y k T s k ) C k ( I − y k s k T y k T s k ) + s k s k T y k T s k C_{k+1}=(I-\frac{s_ks_k^T}{y_k^Ts_k})C_k(I- \frac{y_ks_k^T}{y_k^Ts_k})+\frac{s_ks_k^T}{y_k^Ts_k} Ck+1=(I−ykTskskskT)Ck(I−ykTskykskT)+ykTskskskT(8.29)
存储海森矩阵(Hessian)需要消耗 O ( D 2 ) O(D^2) O(D2)的存储空间,所有对于很大规模的问题,可以使用限制内存BFGS算法(limited memory BFGS),缩写为L-BFGS,其中的 H k H_k Hk或者 H k − 1 H_k^{-1} Hk−1都是用对角矩阵加上低阶矩阵来近似的.具体来说就是积(product) H k − 1 g k H_k^{-1}g_k Hk−1gk可以通过一系列的 s k s_k sk和 y k y_k yk的内积来得到,只使用m个最近的 s k , y k s_k,y_k sk,yk对,忽略掉更早的信息.这样存储上就只需要 O ( m D ) O(mD) O(mD)规模的空间了.通常设置m大约在20左右,就足够有很好的性能表现了.更多相关信息参考(Nocedal and Wright 2006, p177).L-BFGS在机器学习领域中的无约束光滑优化问题中通常都是首选方法.
8.3.6 l 2 l_2 l2规范化(regularization)
相比之下我们更倾向于选岭回归而不是线性回归,类似得到了,对于逻辑回归我们更应该选最大后验估计(MAP)而不是计算最大似然估计(MLE).实际上即便数据规模很大了,在分类背景下规范化还是很重要的.假设数据是线性可分的.这时候最大似然估计(MLE)就可以通过 ∣ ∣ m ∣ ∣ → ∞ ||m||\rightarrow \infty ∣∣m∣∣→∞来得到,对应的就是一个无穷陡峭的S形函数(sigmoid function) I ( w T x > w 0 ) I(w^Tx>w_0) I(wTx>w0),这也叫一个线性阈值单元(linear threshold unit).这将训练数据集赋予了最大规模概率质量.不过这样一来求解就很脆弱而且不好泛化.
所以就得和岭回归里面一样使用 l 2 l_2 l2规范化(regularization).这样一来新的目标函数.梯度函数.海森矩阵就如下所示:
f ′ ( w ) = N L L ( w ) + λ w T w (8.30) g ′ ( w ) = g ( w ) + λ w (8.31) H ′ ( w ) = H ( w ) + λ I (8.32) \begin{aligned} f'(w)&=NLL(w)+\lambda w^Tw &\text{(8.30)}\\ g'(w)&= g(w)+\lambda w &\text{(8.31)}\\ H'(w)&= H(w)+\lambda I &\text{(8.32)}\\ \end{aligned} f′(w)g′(w)H′(w)=NLL(w)+λwTw=g(w)+λw=H(w)+λI(8.30)(8.31)(8.32)
调整过之后的等式就更适合用于各种基于梯度的优化器了.
8.3.7 多类逻辑回归
接着就要讲多类逻辑回归(multinomial logistic regression)了,也叫作最大熵分类器(maximum entropy classifier).这种方法用到的模型形式为:
p ( y = c ∣ x , W ) = exp ( w c T x ) ∑ c ′ = 1 C exp ( w c ′ T x ) p(y=c|x,W)=\frac{\exp(w_c^Tx)}{\sum^C_{c'=1}\exp(w_{c'}^Tx)} p(y=c∣x,W)=∑c′=1Cexp(wc′Tx)exp(wcTx)(8.33)
还有一个轻微变体,叫做条件逻辑模型(conditional logit model),是在对每个数据情况下一系列不同的类集合上进行了规范化;这个可以用于对用户在不同组合的项目集合之间进行的选择进行建模.
然后引入记号.设 μ i c = p ( y i = c ∣ x i W ) = S ( η i ) c \mu_{ic}=p(y_i=c|x_iW)=S(\eta_i)c μic=p(yi=c∣xiW)=S(ηi)c,其中的 η i W T x i \eta_iW^Tx_i ηiWTxi是一个 C × 1 C\times 1 C×1 向量.然后设 y i c = I ( y i = c ) y_{ic}=I(y_i=c) yic=I(yi=c)是对 y i y_i yi的一种编码方式(one-of-C encoding);这样使得 y i y_i yi成为二进制位向量(bit vector),当且仅当 y i = c y_i=c yi=c的时候第c个位值为1.接下来是参考(Krishnapuram et al. 2005),设 w C = 0 w_C=0 wC=0,这样能保证可辨识性(identifiability),然后定义 w = v e c ( W ( : , 1 : C − 1 ) ) w=vec(W(:,1:C-1)) w=vec(W(:,1:C−1))是一个 D × ( C − 1 ) D\times (C-1) D×(C−1)维度的列向量(column vector).
这样就可以写出如下面形式的对数似然函数(log-likelihood):
l ( W ) = log ∏ i = 1 N ∏ c = 1 C μ i c y i c = ∑ i = 1 N ∑ c = 1 C y i c log μ i c (8.34) = ∑ i = 1 N [ ( ∑ c = 1 C y i c w c T x i ) − log ( ∑ c ′ = 1 C exp ( w c ′ T x i ) ) ] (8.35) \begin{aligned} l(W)&= \log\prod^N_{i=1}\prod^C_{c=1}\mu_{ic}^{y_{ic}}=\sum^N_{i=1}\sum^C_{c=1}y_{ic}\log \mu_{ic} &\text{(8.34)}\\ &= \sum^N_{i=1}[(\sum^C_{c=1}y_{ic}w_c^Tx_i)-\log(\sum^C_{c'=1} \exp (w_{c'}^Tx_i) )] &\text{(8.35)}\\ \end{aligned} l(W)=logi=1∏Nc=1∏Cμicyic=i=1∑Nc=1∑Cyiclogμic=i=1∑N[(c=1∑CyicwcTxi)−log(c′=1∑Cexp(wc′Txi))](8.34)(8.35)
加上负号就是负对数似然函数NLL了:
f
(
w
)
=
−
l
(
w
)
f(w)=-l(w)
f(w)=−l(w)
然后就可以针对NLL计算器梯度和海森矩阵了.由于w是分块结构的(block-structured),记号会比较麻烦,但思路还是很简单的.可以定义一个 A ⊗ B A\otimes B A⊗B表示在矩阵A和B之间的克罗内克积(kronecker product).如果A是一个 m × n m\times n m×n矩阵,B是一个 p × q p\times q p×q矩阵,那么 A ⊗ B A\otimes B A⊗B就是一个 m p × n q mp\times nq mp×nq矩阵:
A ⊗ B = [ a 11 B . . . a 1 n B . . . . . . . . . a m 1 B . . . a m n B ] A\otimes B= \begin{bmatrix} a_{11}B &...& a_{1n}B\\...&...&...\\a_{m1}B&...&a_{mn}B \end{bmatrix} A⊗B= a11B...am1B.........a1nB...amnB (8.37)
回到咱们刚刚的例子,很明显(练习8.4)梯度为:
g
(
W
)
=
∇
f
(
w
)
=
∑
i
=
1
N
(
μ
i
−
y
i
)
⊗
x
i
g(W)=\nabla f(w)=\sum^N_{i=1}(\mu_i-y_i)\otimes x_i
g(W)=∇f(w)=∑i=1N(μi−yi)⊗xi(8.38)
其中的 y i = ( I ( y i = 1 ) , . . . , I ( y i = C − 1 ) ) y_i=(I(y_i=1),...,I(y_i=C-1)) yi=(I(yi=1),...,I(yi=C−1)), μ i ( W ) = [ p ( y i = 1 ∣ x i , W ) , . . . , p ( y i = C − 1 ∣ x i , W ) ] \mu_i(W)=[p(y_i=1|x_i,W),...,p(y_i=C-1|x_i,W)] μi(W)=[p(yi=1∣xi,W),...,p(yi=C−1∣xi,W)]这两个都是长度为 C − 1 C-1 C−1的列向量(column vectors).例如如果有特征维度D=3,类别数目C=3,这就成了:
g ( W ) = ∑ i ( ( μ i 1 − y i 1 ) x i 1 ( μ i 1 − y i 1 ) x i 2 ( μ i 1 − y i 1 ) x i 3 ( μ i 2 − y i 2 ) x i 1 ( μ i 2 − y i 2 ) x i 2 ( μ i 2 − y i 2 ) x i 3 ) g(W)=\sum_i\begin{pmatrix}(\mu_{i1}-y_{i1})x_{i1} \\(\mu_{i1}-y_{i1})x_{i2}\\(\mu_{i1}-y_{i1})x_{i3}\\(\mu_{i2}-y_{i2})x_{i1}\\(\mu_{i2}-y_{i2})x_{i2}\\(\mu_{i2}-y_{i2})x_{i3}\end{pmatrix} g(W)=∑i (μi1−yi1)xi1(μi1−yi1)xi2(μi1−yi1)xi3(μi2−yi2)xi1(μi2−yi2)xi2(μi2−yi2)xi3 (8.39)
也就是说对于每个类c,第c列中权重的导数(the derivative for the weights)为:
∇
w
c
f
(
W
)
=
∑
i
(
μ
i
i
−
y
i
c
)
x
i
\nabla_{w_c}f(W)=\sum_i(\mu_i{i}-y_{ic})x_i
∇wcf(W)=∑i(μii−yic)xi(8.40)
这就和二值化逻辑回归的情况下形式一样了,名义上就是误差项乘以 x i x_i xi.(其实这是指数族分布的一个通用特征,在本书9.3.2会讲到.)
(参考练习8.4)很容易发现海森矩阵是一个下面所示形式的 D ( C − 1 ) × D ( C − 1 ) D(C-1)\times D(C-1) D(C−1)×D(C−1)分块矩阵:
$H(WW) = \nabla^2 f(w)=\sumN_{i=1}(diag(\mu_i)-\mu_i\mu_iT)\otimes (x_ix_i^T) $(8.41)
H ( W ) = ∑ i ( μ i 1 − μ i 1 2 − μ i 1 μ i 2 − μ i 1 μ i 2 μ i 2 − μ i 2 2 ) ⊗ ( x i 1 x i 1 x i 1 x i 2 x i 1 x i 3 x i 2 x i 1 x i 2 x i 2 x i 2 x i 3 x i 3 x i 1 x i 3 x i 2 x i 3 x i 3 ) (8.42) = ∑ i ( ( μ i 1 − μ i 1 2 ) X i − μ i 1 μ i 2 X i − μ i 1 μ i 2 X i ( μ i 2 − μ i 2 2 ) X i ) (8.43) \begin{aligned} H(W)&=\sum_i \begin{pmatrix} \mu_{i1}-\mu_{i1}^2 & -\mu_{i1}\mu_{i2}\\ -\mu_{i1}\mu_{i2}& \mu_{i2}-\mu_{i2}^2\end {pmatrix} \otimes \begin{pmatrix} x_{i1}x_{i1}&x_{i1}x_{i2}&x_{i1}x_{i3}\\x_{i2}x_{i1}&x_{i2}x_{i2}&x_{i2}x_{i3}\\x_{i3}x_{i1}&x_{i3}x_{i2}&x_{i3}x_{i3} \end{pmatrix} &\text{(8.42)}\\ &= \sum_i\begin{pmatrix} (\mu_{i1}-\mu_{i1}^2)X_i&-\mu_{i1}\mu_{i2}X_i\\-\mu_{i1}\mu_{i2}X_i&(\mu_{i2}-\mu_{i2}^2)X_i \end{pmatrix} &\text{(8.43)}\\ \end{aligned} H(W)=i∑(μi1−μi12−μi1μi2−μi1μi2μi2−μi22)⊗ xi1xi1xi2xi1xi3xi1xi1xi2xi2xi2xi3xi2xi1xi3xi2xi3xi3xi3 =i∑((μi1−μi12)Xi−μi1μi2Xi−μi1μi2Xi(μi2−μi22)Xi)(8.42)(8.43)
其中的 X i = x i x i T X_i=x_ix_i^T Xi=xixiT.也就是分块矩阵中块 c , c ′ c,c' c,c′部分为:
H c , c ′ ( W ) = ∑ i μ i c ( δ c , c ′ − μ i , c ′ ) x i x i T H_{c,c'}(W)=\sum_i\mu_{ic}(\delta_{c,c'}-\mu_{i,c'})x_ix_i^T Hc,c′(W)=∑iμic(δc,c′−μi,c′)xixiT(8.44)
这也是一个正定矩阵,所以有唯一最大似然估计(MLE).
接下来考虑最小化下面这个式子:
f
′
(
w
W
=
△
−
log
[
(
D
∣
w
)
−
log
p
(
W
)
f'(wW\overset{\triangle}{=} -\log [(D|w)-\log p(W)
f′(wW=△−log[(D∣w)−logp(W)(8.45)
其中的
p
(
W
)
=
∏
c
N
(
w
C
∣
0
,
V
0
)
p(W)=\prod_cN(w_C|0,V_0)
p(W)=∏cN(wC∣0,V0).这样一来就得到了下面所示的新的目标函数/梯度函数/海森矩阵:
f
′
(
w
)
=
f
(
w
)
+
1
2
∑
c
w
c
V
0
−
1
w
c
(8.46)
g
′
(
w
)
=
g
(
W
)
+
V
0
−
1
(
∑
c
w
c
)
(8.47)
H
′
(
w
)
=
H
(
W
)
+
I
C
⊗
V
0
−
1
(8.48)
\begin{aligned} f'(w)&= f(w)+\frac{1}{2}\sum_cw_cV_0^{-1}w_c &\text{(8.46)}\\ g'(w)&= g(W)+V_0^{-1}(\sum_cw_c) &\text{(8.47)}\\ H'(w)&= H(W)+I_C\otimes V_0^{-1} &\text{(8.48)}\\ \end{aligned}
f′(w)g′(w)H′(w)=f(w)+21c∑wcV0−1wc=g(W)+V0−1(c∑wc)=H(W)+IC⊗V0−1(8.46)(8.47)(8.48)
这样就可以传递给任意的基于梯度的优化器来找到最大后验估计(MAP)了.不过要注意这时候的海森矩阵规模是 O ( ( C D ) × ( C D ) ) O ((CD)\times (CD)) O((CD)×(CD)),比二值化情况下要多C倍的行和列,所以这时候更适合使用限制内存的L-BFGS方法,而不适合用牛顿法.具体MATLAB代码参考本书配套PMTK3的logregFit.
8.4 贝叶斯逻辑回归
对逻辑回归模型,很自然的需求就是计算在参数上的完整后验分布
p
(
w
∣
D
)
p(w|D)
p(w∣D).只要我们相对预测指定一个置信区间,这就很有用了.(另外在解决情景匪徒问题(contextual bandit problems)时候也很有用,参考本书5.7.3.1.)
然而很不幸,不像线性回归的时候了,这时候没办法实现这个目的,因为对于逻辑回归来说没有合适的共轭先验.本章这里讲的是一个简单近似;还有一些其他方法,比如马尔可夫链蒙特卡罗(Markov Chain Monte Carlo,缩写为MCMC,本书24.3.3.1),变分推导(variational inference,本书21.8.1.1),
期望传播(expectation propagation,Kuss and Rasmussen 2005)等等.为了记号简单,这里还用二值逻辑回归为例.
8.4.1 拉普拉斯近似(Laplace approximation)
本节将对一个后验分布进行高斯近似.假如
θ
∈
R
D
\theta\in R^D
θ∈RD.设:
p
(
θ
∣
D
)
=
1
Z
e
−
E
(
θ
)
p(\theta|D)=\frac{1}{Z}e^{-E(\theta)}
p(θ∣D)=Z1e−E(θ)(8.49)
其中的
E
(
θ
)
E(\theta)
E(θ)叫能量函数(energy function),等于未归一化对数后验(unnormalized
log posterior)的负对数,即
E
(
θ
)
=
−
log
p
(
θ
,
D
)
E(\theta)=-\log p(\theta,D)
E(θ)=−logp(θ,D),
Z
=
p
(
D
)
Z=p(D)
Z=p(D)是归一化常数.然后进行关于众数
θ
∗
\theta^*
θ∗(对应最低能量状态)的泰勒级数展开,就得到了:
E ( θ ) ≈ E ( θ ∗ ) + ( θ − θ ∗ ) T g + 1 2 ( θ − θ ∗ ) T H ( θ − θ ∗ ) E(\theta) \approx E(\theta^*)+(\theta-\theta^*)^Tg+\frac{1}{2}(\theta-\theta^*)^TH(\theta-\theta^*) E(θ)≈E(θ∗)+(θ−θ∗)Tg+21(θ−θ∗)TH(θ−θ∗)(8.50)
其中的g是梯度,H是在众数位置能量函数的海森矩阵:
g
=
△
∇
E
(
θ
)
∣
θ
∗
,
H
=
△
∂
2
E
(
θ
)
∂
θ
∂
θ
T
∣
θ
∗
g\overset{\triangle}{=} \nabla E(\theta)|_{\theta^*},H\overset{\triangle}{=} \frac{\partial^2E(\theta)}{\partial\theta\partial\theta^T}|_{\theta^*}
g=△∇E(θ)∣θ∗,H=△∂θ∂θT∂2E(θ)∣θ∗(8.51)
由于
θ
∗
\theta^*
θ∗是众数,梯度项为零.因此:
p
^
(
θ
∣
D
)
≈
1
Z
e
−
E
(
θ
∗
)
exp
[
−
1
2
(
θ
−
θ
∗
)
T
H
(
θ
−
θ
∗
)
]
(8.52)
=
N
(
θ
∣
θ
∗
,
H
−
1
)
(8.53)
Z
=
p
(
D
)
≈
∫
p
^
(
θ
∣
D
)
d
θ
=
e
−
E
(
θ
∗
)
(
2
π
)
D
/
2
∣
H
∣
−
1
2
(8.54)
\begin{aligned} \hat p(\theta|D)& \approx \frac{1}{Z}e^{-E(\theta^*)} \exp[-\frac{1}{2}(\theta-\theta^*)^T H(\theta-\theta^*)] &\text{(8.52)}\\ & = N(\theta|\theta^*,H^{-1})&\text{(8.53)}\\ Z=p(D)& \approx \int \hat p(\theta|D)d\theta = e^{-E(\theta^*)}(2\pi)^{D/2}|H|^{-\frac{1}{2}} &\text{(8.54)}\\ \end{aligned}
p^(θ∣D)Z=p(D)≈Z1e−E(θ∗)exp[−21(θ−θ∗)TH(θ−θ∗)]=N(θ∣θ∗,H−1)≈∫p^(θ∣D)dθ=e−E(θ∗)(2π)D/2∣H∣−21(8.52)(8.53)(8.54)
最后这一行是参考了多元高斯分布的归一化常数.
等式8.54就是对边缘似然函数的拉普拉斯近似.所以等式8.52有时候也叫做对后验的拉普拉斯近似.不过在统计学领域,拉普拉斯近似更多指的是一种复杂方法(具体细节参考Rue等,2009).高斯近似通常就足够近似了,因为后验分布随着样本规模增长就越来越高斯化,这个类似中心极限定理.(在物理学领域有一个类似的技术叫做鞍点近似(saddle point approximation).)
8.4.2 贝叶斯信息量(Bayesian information criterio,缩写为BIC)的推导
可以使用高斯近似来写出对数似然函数,去掉不相关的常数之后如下所示:
log
p
(
D
)
≈
log
(
p
(
D
∣
θ
∗
)
+
log
p
(
θ
∗
)
−
1
2
log
∣
H
∣
\log p(D)\approx \log(p(D|\theta^*)+\log p(\theta^*)-\frac{1}{2}\log|H|
logp(D)≈log(p(D∣θ∗)+logp(θ∗)−21log∣H∣(8.55)
加在 log ( p ( D ∣ θ ∗ ) \log(p(D|\theta^*) log(p(D∣θ∗)后面的惩罚项(penalization terms)也叫作奥卡姆因子(Occam factor),是对模型复杂程度的量度.如果使用均匀先验,即 p ( θ ) ∝ 1 p(\theta)\propto 1 p(θ)∝1,就可以去掉第二项,然后把 θ ∗ \theta^* θ∗替换成最大似然估计(MLE) θ ^ \hat\theta θ^.
接下来看看对上式中第三项的估计.已知 H = ∑ i = 1 N H i H=\sum^N_{i=1}H_i H=∑i=1NHi,其中的 H i = ∇ ∇ log p ( D i ∣ θ ) H_i=\nabla\nabla \log p(D_i|\theta) Hi=∇∇logp(Di∣θ).然后用一个固定的矩阵 H ^ \hat H H^来近似每个 H i H_i Hi.这样就得到了:
log ∣ H ∣ = log ∣ N H ^ ∣ = log ( N d ∣ H ^ ∣ ) = D log N + log ∣ H ^ ∣ \log|H|=\log|N\hat H| =\log (N^d|\hat H|)=D\log N+\log |\hat H| log∣H∣=log∣NH^∣=log(Nd∣H^∣)=DlogN+log∣H^∣(8.56)
其中
D
=
d
i
m
(
θ
)
D=dim(\theta)
D=dim(θ),并且假设H是满秩矩阵(full rank matrix).就可以去掉
log
∣
H
^
∣
\log |\hat H|
log∣H^∣这一项了,因为这个独立于N,这样就可以被似然函数盖过去了.将所有条件结合起来,就得到了贝叶斯信息量分数(BIC score,参考本书5.3.2.4):
log
p
(
D
)
≈
p
(
D
∣
θ
^
)
−
D
2
log
N
\log p(D)\approx p(D|\hat \theta)-\frac{D}{2}\log N
logp(D)≈p(D∣θ^)−2DlogN(8.57)
8.4.3 逻辑回归的高斯近似
接下来对逻辑回归应用高斯近似.使用一个高斯先验,形式为
p
(
w
)
=
N
(
w
∣
0
,
V
0
)
p(w)=N(w|0,V_0)
p(w)=N(w∣0,V0),正如在最大后验估计(MAP)里面一样.近似后验为:
p
(
w
∣
D
)
≈
N
(
w
∣
w
^
,
H
−
1
)
p(w|D)\approx N(w|\hat w, H^{-1})
p(w∣D)≈N(w∣w^,H−1)(8.58)
上式中 w ^ = arg min w E ( w ) , E ( w ) = − ( log p ( D ∣ w ) + log p ( w ) ) , H = ∇ 2 E ( w ) ∣ w ^ \hat w = \arg\min_w \mathrm{E}(w),\mathrm{E}(w) =-(\log p(D|w)+\log p(w)), H=\nabla^2 \mathrm{E}(w)|_{\hat w} w^=argminwE(w),E(w)=−(logp(D∣w)+logp(w)),H=∇2E(w)∣w^
以图8.5(a)中所示的线性可分二位数据为例.有很多个不同参数设置的线都能很好地区分开训练数据;在图中画了四条.图8.5(b)所示的是似然函数曲面,从中可以看到似然函数是无界的,随着向右上角的参数空间移动,似然函数沿着一个嵴线移动, w 2 / w 1 = 2.34 w_2/w_1=2.34 w2/w1=2.34(通过对角线推测到的).这是因为将 ∣ ∣ w ∣ ∣ ||w|| ∣∣w∣∣推向无穷大可以是似然函数最大化.因为大的回归权重会让S形函数非常陡峭,就成了一个阶梯函数了.因此当数据线性可分的时候最大似然估计(MLE)不能很好定义.
为了规范表达这个问题,假设使用一个以原点为中心的模糊球形先验(vague spherical prior), N ( w ∣ 0 , 100 I ) N(w|0,100I) N(w∣0,100I).将这个球形先验乘以似然函数曲面就得到了一个严重偏斜的后验,如图8.5©所示.(这是因为似然函数截断了参数空间中与实验数据不符合的区域.)最大后验估计(MAP)如图中蓝色点所示.这就不像是最大似然估计(MLE)那样跑到无穷远了.
对这个后验的高斯近似如图8.5(d)所示.其中可以看到是一个对称分布,因此也不算是一个很好的估计.不过好歹众数还是正确的(通过构造得到的),而且至少表现了沿着西南东北方向比垂直这个方向上有更大不确定性的这个事实(这对应着稀疏线方向的不确定性).虽然这个高斯近似很粗糙,也肯定比用最大后验估计(MAP)得到的 δ \delta δ函数近似要好很多了.
8.4.4 近似后验预测
给定了后验,就可以计算置信区间,进行假设检验等等.就跟之前在7.6.3.3里面对线性回归案例所做的一样.不过在机器学习里面,更多兴趣还是在预测上面.后验预测分布的形式为:
p
(
y
∣
x
,
D
)
=
∫
p
(
y
∣
x
,
w
)
p
(
w
∣
D
)
d
w
p(y|x,D)=\int p(y|x,w)p(w|D)dw
p(y∣x,D)=∫p(y∣x,w)p(w∣D)dw(8.59)
此处参考原书图8.5
很不幸这个积分可难算了.
最简单的近似就是插值估计(plug-in approximation),在二值化分类的情况下,形式为:
p
(
y
=
1
∣
x
,
D
)
≈
p
(
y
=
1
∣
x
,
E
[
w
]
)
p(y=1|x,D)\approx p(y=1|x,\mathrm{E}[w])
p(y=1∣x,D)≈p(y=1∣x,E[w])(8.60)
其中的 E [ w ] \mathrm{E}[w] E[w]是后验均值.在这个语境下, E [ w ] \mathrm{E}[w] E[w]也叫作贝叶斯点(Bayes point).当然了,这种插值估计肯定是低估了不确定性了.所以接下来会说更好的近似方法.
此处参考原书图8.6
8.4.4.1 蒙特卡罗方法近似
更好的方法就是蒙特卡罗方法(Monte Carlo approximation),定义如下:
p
(
y
=
1
∣
x
,
D
)
≈
1
S
∑
s
=
1
S
s
i
g
m
(
(
w
s
)
T
x
)
p(y=1|x,D)\approx \frac{1}{S}\sum^S_{s=1} sigm((w^s)^Tx)
p(y=1∣x,D)≈S1∑s=1Ssigm((ws)Tx)(8.61)
其中的 w s ∼ p ( w ∣ D ) w^s\sim p(w|D) ws∼p(w∣D)是在后验中的取样.(这个方法很容易扩展到多类情况.)如果使用蒙特卡罗方法估计后验,就可以服用这些样本来进行预测.如果对后验使用高斯估计,就要用标准方法从高斯分布中取得独立样本.
此处查看原书图8.7
图8.6(b)展示了在我们的二维样例中从后验预测分布取的样本.图8.6©展示的是这些样本的均值.通过对多个预测取平均值,就可以发现决策边界的不确定性随着远离训练数据而散发开.所以虽然决策边界是线性的,后验预测密度还是非线性的.要注意后验均值的决策边界大概和每个类的距离都是一样远;这是对大便捷原则的贝叶斯模拟,在本书14.5.2.2会进一步讲到.
图8.7(a)所示的是一维下第一个例子.红色的点表示的是训练数据中后验预测分布评估出来的均值.竖着的蓝色线段表示的是对该后验预测的95%置信区间;蓝色的小星星是中位数位置.可见贝叶斯方法可以基于SAT分数来对一个学生能通过考试的概率的不确定性进行建模,而不是仅仅给出一个点估计.
8.4.4.2 Probit近似(适度输出)*
如果有一个对后验的高斯近似, p ( w ∣ D ) ≈ N ( w ∣ m N , V N ) p(w|D) \approx N(w|m_N,V_N) p(w∣D)≈N(w∣mN,VN),就可以计算后验预测分布的确定性近似(deterministic approximation),至少在二值化情况下是可以的.步骤如下所示:
p ( y = 1 ∣ x , D ) ≈ ∫ s i g m ( w T x ) p ( w ∣ D ) d w = ∫ s i g m a ( a ) N ( a ∣ μ a , σ a 2 ) d a (8.62) a = △ w T x (8.63) μ a = △ E [ a ] = m N T x (8.64) σ a 2 = △ v a r [ a ] = ∫ p ( a ∣ D ) [ a 2 − E [ a 2 ] ] d a (8.65) = ∫ p ( w ∣ D ) [ ( w T x ) 2 − ( m N T x ) 2 ] d w = x T V N x (8.66) \begin{aligned} p(y=1|x,D)&\approx \int sigm(w^Tx)p(w|D)dw=\int sigma(a)N(a|\mu_a,\sigma_a^2)da &\text{(8.62)}\\ a&\overset{\triangle}{=} w^Tx &\text{(8.63)}\\ \mu_a &\overset{\triangle}{=} \mathrm{E}[a]=m_N^Tx &\text{(8.64)}\\ \sigma^2_a&\overset{\triangle}{=} var[a]=\int p(a|D)[a^2-\mathrm{E}[a^2]]da &\text{(8.65)}\\ &= \int p(w|D)[(w^Tx)^2-(m_N^Tx)^2]dw =x^TV_Nx &\text{(8.66)}\\ \end{aligned} p(y=1∣x,D)aμaσa2≈∫sigm(wTx)p(w∣D)dw=∫sigma(a)N(a∣μa,σa2)da=△wTx=△E[a]=mNTx=△var[a]=∫p(a∣D)[a2−E[a2]]da=∫p(w∣D)[(wTx)2−(mNTx)2]dw=xTVNx(8.62)(8.63)(8.64)(8.65)(8.66)
然后就能发现需要去评估一个对应高斯分布的S形函数(sigmoid function)的期望.可以利用S形函数类似概率函数(probit function)的性质来进行近似,通过标准正态分布的累计密度函数(cdf):
Φ
(
a
)
=
△
∫
−
∞
a
N
(
x
∣
0
,
1
)
d
x
\Phi(a)\overset{\triangle}{=}\int^a_{-\infty} N(x|0,1)dx
Φ(a)=△∫−∞aN(x∣0,1)dx(8.67)
图8.7(b)所示的就是s形函数和概率函数.图中的坐标轴做了缩放,使得 s i g m ( a ) sigm(a) sigm(a)在原点位置附近有和 Φ ( λ a ) \Phi(\lambda a) Φ(λa)有类似的范围,其中 λ 2 = π / 8 \lambda ^2 =\pi /8 λ2=π/8.
利用概率函数(probit function)的优势就是可以用一个高斯分布以解析形式卷积(convolve)出来:
∫
Φ
(
λ
a
)
N
(
a
∣
μ
,
σ
2
)
d
a
=
Φ
(
a
(
λ
−
2
+
σ
2
)
1
2
)
\int \Phi (\lambda a)N(a|\mu,\sigma^2)da =\Phi (\frac{a}{(\lambda ^{-2}+\sigma^2)^{\frac{1}{2}}})
∫Φ(λa)N(a∣μ,σ2)da=Φ((λ−2+σ2)21a)(8.68)
然后再等号两边都使用 s i g m ( a ) ≈ Φ ( λ a ) sigm(a)\approx \Phi (\lambda a) sigm(a)≈Φ(λa) 来插值近似,就得到了:
∫ s i g m ( a ) N ( a ∣ μ , σ 2 ) d a ≈ s i g m ( k ( σ 2 ) μ ) (8.69) k ( σ 2 ) = △ ( 1 + π σ 2 / 8 ) − 1 2 (8.70) \begin{aligned} \int sigm(a)N(a|\mu,\sigma^2)da& \approx sigm(k(\sigma^2)\mu) &\text{(8.69)}\\ k(\sigma^2)& \overset{\triangle}{=} (1+\pi\sigma^2/8)^{-\frac{1}{2}} &\text{(8.70)}\\ \end{aligned} ∫sigm(a)N(a∣μ,σ2)dak(σ2)≈sigm(k(σ2)μ)=△(1+πσ2/8)−21(8.69)(8.70)
将上面的近似用到逻辑回归模型中,就得到了下面的表达式(最初引用自(Spiegelhalter and Lauritzen 1990)):
p
(
y
=
1
∣
x
,
D
)
≈
s
i
g
m
(
k
(
σ
a
2
)
μ
a
)
p(y=1|x,D)\approx sigm(k(\sigma^2_a)\mu_a)
p(y=1∣x,D)≈sigm(k(σa2)μa)(8.71)
图8.6(d)所示表面这个近似能得到和蒙特卡洛近似相似的结果.
使用等式8.71这种近似,也叫作调节输出(moderated output),因为这比起插值估计不那么极端.要理解这一点,要注意到
0
≤
k
(
σ
2
)
≤
1
0\le k(\sigma^2)\le 1
0≤k(σ2)≤1,因此就有:
s
i
g
m
(
k
(
σ
2
)
μ
)
≤
s
i
g
m
a
(
μ
)
=
p
(
y
=
1
∣
x
,
w
^
)
sigm(k(\sigma^2)\mu)\le sigma(\mu)=p(y=1|x,\hat w)
sigm(k(σ2)μ)≤sigma(μ)=p(y=1∣x,w^)(8.72)
上式中的不等号当
μ
≠
0
\mu\ne 0
μ=0的时候严格成立.如果
μ
>
0
\mu >0
μ>0,就有
p
(
y
=
1
∣
x
,
w
^
)
>
0.5
p(y=1|x,\hat w)>0.5
p(y=1∣x,w^)>0.5,但调节的预测(moderated prediction)总是在0.5附近,所以不太可信(less confident).不过,决策边界存在于
p
(
y
=
1
∣
x
,
D
)
=
s
i
g
m
(
k
(
σ
2
)
μ
)
=
0.5
p(y=1|x,D)=sigm(k(\sigma^2)\mu)=0.5
p(y=1∣x,D)=sigm(k(σ2)μ)=0.5的位置,也就意味着
μ
=
w
^
T
x
=
0
\mu =\hat w^Tx=0
μ=w^Tx=0.因此调节预测的决策边界和差值近似是一样的.所以两种方法的误分类率是一样的,但对数似然函数是不一样的.(要注意在多类情况下,后验协方差给出的结果和插值方法是不一样的,参考本书练习3.10.3以及(Rasmussen
and Williams 2006)).
8.4.5 残差分析(异常值检测)*
检验数据中的异常值(outlier)有时候很有用.这个过程就叫残差分析(residual analysis)或者案例分析(case analysis).在回归情况下,这个可以通过计算残差来得到: r i = y i − y ^ i r_i=y_i-\hat y_i ri=yi−y^i,其中的 y ^ i = w ^ T x i \hat y_i =\hat w^T x_i y^i=w^Txi.如果模型假设正确无误,这些值应该遵循一个正态分布 N ( 0 , σ 2 ) N(0,\sigma^2) N(0,σ2).可以投图(qq-plot)评定,其中两个坐标轴的值分别是高斯分布的N个理论值(theoretical quantiles)和 r i r_i ri的经验值(empirical quantiles).偏离这条直线的点就是潜在的异常值.
基于残差的简单方法不适用于二进制数据,因为依靠测试统计的渐进正态(asymptotic normality).不过用贝叶斯方法,就能定义异常值了,可以定义为 p ( y i ∣ y ^ i ) p(y_i |\hat y_i) p(yi∣y^i)小的点,一般使用 y ^ i = s i g m ( w ^ T x ) \hat y_i =sigm (\hat w^T x) y^i=sigm(w^Tx).要注意这里的 w ^ \hat w w^是从全部数据估计出来的.在预测 y i y_i yi时更好的方法是从对w的估计排除掉 ( x i , y i ) (x_i,y_i) (xi,yi).也就是定义异常值为在交叉验证后验预测分布下有低概率的点,定义形式为:
p ( y i ∣ x i , x − i , y − i ) = ∫ p ( y i ∣ x i , w ) ∏ i ′ ≠ i p ( y i ′ ∣ x i ′ , w ) p ( w ) d w p(y_i|x_i,x_{-i},y_{-i})=\int p(y_i|x_i,w)\prod _{i'\ne i} p(y_{i'}|x_{i'},w)p(w)dw p(yi∣xi,x−i,y−i)=∫p(yi∣xi,w)∏i′=ip(yi′∣xi′,w)p(w)dw(8.73)
这就可以通过抽样方法来有效近似(Gelfand 1996).更多关于逻辑回归模型中残差分析的相关内容可以参考(Johnson and Albert 1999, Sec 3.4).
扫一扫
6768
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
