09 Generalized linear models and the exponential family

9.3 广义线性模型(Generalized linear models,缩写为GLMs)

线性回归和逻辑回归都属于广义线性模型的特例(McCullagh and Nelder 1989).
这些模型中输出密度都是指数族分布(参考本书9.2),而均值参数都是输入的线性组合,经过可能是非线性的函数,比如逻辑函数等等.下面就要详细讲一下广义线性模型(GLMs).为了记号简单,先看标量输出的情况.(这就排除了多远逻辑回归,不过这只是为表述简单而已.)

9.3.1 基础知识

要理解广义线性模型,首先要考虑一个标量响应变量的无条件分布(unconditional dstribution)的情况:

$p(y_i|\theta ,{\sigma }^2)=\exp [\frac{y_i\theta -A(\theta )}{{\sigma }^2}+c(y_i,{\sigma }^2)]$(9.77)

上式中的${\sigma }^2$叫做色散参数(dispersion parameter),通常设为1.$\theta$是自然参数,A是配分函数,c是归一化常数.例如,在逻辑回归的情况下,$\theta$就是对数比值比(log-odds ratio),$\theta =\log (\frac{\mu }{1-\mu })$,其中$\mu =\mathrm{E}[y]=p(y=1)$是均值参数(mean parameter),参考本书9.2.2.1.要从均值参数转成自然参数(natural parameter),可以使用一个函数$\varphi$,也就是$\theta =\Psi (\mu )$.这个函数由指数族分布的形式唯一确定(uniquely determined).实际上这是一个可逆映射(invertible mapping),所以也就有$\mu ={\Psi }^{-1}(\theta )$.另外通过本书9.2.3可以知道这个均值可以通过对配分函数(partition function)求导而得到,也就是有$\mu ={\Psi }^{-1}(\theta )=A^{\prime }(\theta )$.

然后加上输入/协变量(covariates).先定义一个输入特征的线性函数:

${\eta }_i=w^T{x}_i$(9.78)

分布	Link $g(\mu )$	$\theta =\psi (\mu )$	$\mu ={\psi }^{-1}(\theta )=\mathrm{E}[y]$
$N(\mu ,{\sigma }^2)$	indentity	$\theta =\mu$	$\mu =\theta$
$Bin(N,\mu )$	logit	$\theta =\log \frac{\mu }{1-\mu }$	$\mu =sigm(\theta )$
$Poi(\mu )$	log	$\theta =\log (\mu )$	$\mu =e^{\theta }$

表 9.1 常见广义线性模型(GLMs)的连接函数(link function)$\psi$.

然后使这个分布的均值为这个线性组合的某个可逆单调函数.通过转换,得到这个函数就叫做均值函数(mean function),记作$g^{-1}$,所以:
${\mu }_i=g^{-1}({\eta }_i)=g^{-1}(w^T{x}_i)$(9.79)

如图9.1所示为这个简单模型的总结.

均值函数(mean function)的逆函数,记作$g()$,就叫做连接函数(link function).我们可以随意选择任意函数来作为连接函数,只要是可逆的,以及均值函数$g()$有适当的范围.例如在逻辑回归里面,就设置${\mu }_i=g^{-1}({\eta }_i)=sigm({\eta }_i)$.

连接函数有一个特别简单的形式,就是$g=\varphi$,这也叫做规范连接函数(canonical link function).这种情况下则有${\theta }_i={\eta }_i=w^T{x}_i$,所以模型就成了:

$p(y_i|x_i,w,{\sigma }^2)=\exp [\frac{y_i{w}^T{x}_i-A(w^T{x}_i)}{{\sigma }^2}+c(y_i,{\sigma }^2)]$(9.80)

表格9.1中所示的是一些分布和规范连接函数.可见伯努利分布或者二项分布的规范连接函数是$g(\mu )=\log (\eta /(1-\eta ))$,而你函数是逻辑函数(logistic function)$\mu =sigm(\eta )$.

基于本书9.2.3的结果,可以得到响应变量的均值和方差:

$$\begin{array}{rlr}\mathrm{E}[y|x_i,w,{\sigma }^2]& ={\mu }_i=A^{\prime }({\theta }_i)& \text{(9.81)}\\ var[y|x_i,w,{\sigma }^2]& ={\sigma }_i^2=A^{\prime \prime }({\theta }_i){\sigma }^2& \text{(9.82)}\end{array}$$

为了记好清楚,接下来就看一些简单样例.

对于线性回归,则有:

$\log p(y_i|x_i,w,{\sigma }^2)=\frac{y_i{\mu }_i-{\mu }_i^2/2}{{\sigma }^2}-\frac{1}{2}(\frac{y_i^2}{{\sigma }^2}+\log (2\pi {\sigma }^2))$(9.83)

其中$y_i\in R,{\theta }_i={\mu }_i=w^T{x}_i$,而$A(\theta )={\theta }^2/2$,所以$\mathrm{E}[y_i]={\mu }_i,var[y_i]={\sigma }^2$.

对于二项回归(binomial regression),则有:

$\log p(y_i|x_i,w)=y_i\log (\frac{{\pi }_i}{1-{\pi }_i}+N_i\log (1-{\pi }_i)+\log \begin{array}{r}{N}_i\\ y_i\end{array}$(9.84)

其中 y i ∈ { 0 , 1 , . . . , N i } , π i = s i g m ( w T x i ) , θ i = log ⁡ ( π i / ( 1 − π i ) ) = w T x i , σ 2 = 1 y_i \in \{0,1,...,N_i\},\pi_i =sigm(w^Tx_i),\theta_i=\log (\pi_i/(1-\pi_i))=w^Tx_i,\sigma^2=1 yi​∈{0,1,...,Ni​},πi​=sigm(wTxi​),θi​=log(πi​/(1−πi​))=wTxi​,σ2=1. $A(\theta )=N_i\log (1+e^{\theta })$,所以$\mathrm{E}[y_i]=N_i{\pi }_i={\mu }_i,var[y_i]=N_i{\pi }_i(1-{\pi }_i)$

对于泊松分布(poisson regression),则有:
$\log p(y_i|x_i,w)=y_i\log {\mu }_i-{\mu }_i-\log (y_i!)$(9.85)

其中 y i ∈ { 0 , 1 , 2 , . . . } , μ i = exp ⁡ ( w T x i ) , θ = log ⁡ ( μ i ) = w T x i , σ 2 = 1 y_i\in \{0,1,2,...\},\mu_i =\exp (w^Tx_i),\theta=\log(\mu_i)=w^Tx_i,\sigma^2=1 yi​∈{0,1,2,...},μi​=exp(wTxi​),θ=log(μi​)=wTxi​,σ2=1.而$A(\theta )=e^{\theta }$,所以$\mathrm{E}[y_i]=var[y_i]={\mu }_i$.泊松回归在生物统计中应用很广,其中的$y_i$可能代表着给定人群或者地点的病患数目,或者高通量测序背景下基因组位置的读数数量,参考(Kuan et al. 2009).

9.3.2 最大似然估计(MLE)和最大后验估计(MAP)

广义线性模型的最重要的一个性质就是可以用和逻辑回归拟合的同样方法来进行拟合.对数似然函数形式如下所示:

$$\begin{array}{rlr}l(w)=\log p(D|w)& =\frac{1}{{\sigma }^2}\sum _{i=1}^N{l}_i& \text{(9.86)}\\ l_i& \stackrel{△}{=}{\theta }_i{y}_i-A({\theta }_i)& \text{(9.87)}\end{array}$$

$$\begin{array}{rlr}\frac{d{l}_i}{d{w}_j}& =\frac{d{l}_i}{d{\theta }_i}\frac{d{\theta }_i}{d{\mu }_i}\frac{d{\mu }_i}{d{\eta }_i}\frac{d{\eta }_i}{d{w}_j}& \text{(9.88)}\\ & =(y_i-A^{\prime }({\theta }_i))\frac{d{\theta }_i}{d{\mu }_i}\frac{}{}{x}_{ij}& \text{(9.89)}\\ & =(y_i-{\mu }_i)\frac{d{\theta }_i}{d{\mu }_i}\frac{}{}{x}_{ij}& \text{(9.90)}\end{array}$$

名称	公式
Logistic	$g^{-1}(\eta )=sigm(\eta )=\frac{e^{\eta }}{1+e^{\eta }}$
Probit	$g^{-1}(\eta )=\Phi (\eta )$
Log-log	$g^{-1}(\eta )=\exp (-\exp (-\eta ))$
Complementary log-log	$g^{-1}(\eta )=1-\exp (-\exp (\eta ))$
表格9.2 二值回归(binary regression)的一些可能的均值函数的总结.

如果使用规范连接函数${\theta }_i={\eta }_i$,就简化成了:
${\nabla }_{w}l(w)=\frac{1}{{\sigma }^2}[{\sum }_{i=1}^N(y_i-{\mu }_i)x_i]$(9.91)

这就成了输入向量以误差为权重的求和.这可以在一个(很复杂的)梯度下降过程中用到,具体参考本书8.5.2.不过为了提高效率,还是应该用一个二阶方法.如果使用一个规范连接函数,海森矩阵(Hessian)就是:
$H=-\frac{1}{{\sigma }^2}{\sum }_{i=1}^N\frac{d{\mu }_i}{d{\theta }_i}{x}_i{x}_i^T=-\frac{1}{{\sigma }^2}{X}^{T}SX$(9.92)

其中$S=diag(\frac{d{\mu }_1}{d{\theta }_1},...,\frac{d{\mu }_N}{d{\theta }_N})$是对角权重矩阵(diagonal weighting matrix).这可以用在一个迭代重加权最小二乘法(Iteratively reweighted least squares,缩写为IRLS)算法内部(参考本书8.3.4).然后可以有下面所示的牛顿更新:

$$\begin{array}{rlr}{w}_{t+1}& =(X^T{S}_{t}X{)}^{-1}{X}^T{S}_t{z}_t& \text{(9.93)}\\ z_t& ={\theta }_t+S_t^{-1}(y-{\mu }_t)& \text{(9.94)}\end{array}$$

其中${\theta }_t=X{w}_t,{\mu }_t=g^{-1}({\eta }_t)$.

如果我们扩展求导(extend derivation)来处理非规范连接函数(non-canonical links),就会发现海森矩阵(Hessian)有另外一项.不过最终期望海森矩阵(expected Hessian)和等式9.92一模一样,利用这个期望海森矩阵(也就是费舍尔信息矩阵(Fisher information matrix))来替代真实的海森矩阵就叫做费舍尔排序方法(Fisher scoring method).

修改上面的过程来进行高斯先验的最大后验估计(MAP estimation)就很简单了:只需要调整目标函数(objective),梯度*gradient)和海森矩阵(Hessian),就像是在本书8.3.6里面对逻辑回归加上$l_2$规范化(regularization)一样.

9.3.3 贝叶斯推导

广义线性模型(GLMs)的贝叶斯推导通常都用马尔科夫链蒙特卡罗方法(Markov chain Monte Carlo,缩写为 MCMC),参考本书24章.可用基于迭代重加权最小二乘法(IRLS-based proposal)来的Metropolis Hastings方法(Gamerman 1997),或者每个全条件(full conditional)使用自适应拒绝采样(adaptive rejection sampling,缩写为ARS)的吉布斯采样(Gibbs sampling)(Dellaportas and Smith 1993).更多信息参考(Dey et al. 2000).另外也可以用高斯近似(Gaussian approximation,本书8.4.1)或者变分推理(variational inference,本书21.8.11).

9.4 概率单位回归(Probit regression)

在(二值化)洛基回归中,用到了形式为$p(y=1|x_i,w)=sigm(w^T{x}_i)$的模型.对任何从区间$[-\infty ,\infty ]$到$[0,1]$进行映射的函数$g^{-1}$,一般都可以写出$p(y=1|x_i,w)=g^{-1}(w^T{x}_i)$.表格9.2所示为若干可能的均值函数.
在本节,要讲的情况是$g^{-1}(\eta )=\Phi (\eta )$,其中的$\Phi (\eta )$是标准正态分布(standard normal)的累积密度函数(cdf).这样就叫概率单位回归(probit regression).概率函数(probit function)和逻辑函数(logistic function)很相似,如图8.7(b)所示.不过这个模型比逻辑回归有一些优势.

9.4.1 使用基于梯度优化的最大似然估计(MLE)和最大后验估计(MAP)

我们可以对概率单位回归使用标准梯度方法来得到最大似然估计(MLE).设 μ i = w T x i , y ~ i ∈ { − 1 , + 1 } \mu_i=w^Tx_i,\tilde y_i\in \{-1,+1\} μi​=wTxi​,y~​i​∈{−1,+1}.然后对于一个特定情况的对数似然函数的梯度就是:
$g_i\stackrel{△}{=}\frac{d}{dw}\log p({\tilde{y}}_i|w^T{x}_i)=\frac{d{\mu }_i}{dw}\frac{d}{d{\mu }_i}\log p({\tilde{y}}_i)|w^T{x}_i)=x_i\frac{{\tilde{y}}_i\varphi ({\mu }_i)}{\Phi ({\tilde{y}}_i{\mu }_i)}$(9.95)

$H_i=\frac{d}{d{w}^2}\log p({\tilde{y}}_i|w^T{x}_i)=-x_i(\frac{\varphi ({\mu }_i{)}^2}{\Phi ({\tilde{y}}_i{\mu }_i)}+\frac{{\tilde{y}}_i{\mu }_i\varphi ({\mu }_i)}{\Phi ({\tilde{y}}_i{\mu }_i)})x_i^T$(9.96)

$p(w)=N(0,V_0)$ ${\sum }_i{g}_i+2V_0^{-1}w$ ${\sum }_i{H}_i+2V_0^{-1}$

9.4.2 潜在变量解释(latent variable interpretation)

$$\begin{array}{rlr}{u}_{0i}& \stackrel{△}{=}{w}_0^T{x}_i+{\delta }_{0i}& \text{(9.97)}\\ u_{1i}& \stackrel{△}{=}{w}_1^T{x}_i+{\delta }_{1i}& \text{(9.98)}\\ y_i& =I(u_{1i}\ge u_{10})& \text{(9.99)}\end{array}$$

随机效用模型(random utility model,缩写为RUM,McFadden 1974; Train 2009).

${ϵ}_i={\delta }_{1i}-{\delta }_{0i}$

$$\begin{array}{rlr}{z}_i& \stackrel{△}{=}{w}^T{x}_i+{ϵ}_i& \text{(9.100)}\\ {ϵ}_i& \sim N(0,1)& \text{(9.101)}\\ y_i=1& =I(z_i\ge 0)& \text{(9.102)}\end{array}$$

根据(Fruhwirth-Schnatter and Fruhwirth 2010),将这叫做差别随机效用模型(difference random utility model,缩写为dRUM).

当边缘化去掉(marginalize out)$z_i$之后,就恢复(recover)了概率单位模型(probit model):

$$\begin{array}{rlr}p(y_i=1|x_i,w)& =\int I(z_i\ge 0)N(z_i|w^T{x}_i,1)d{z}_i& \text{(9.103)}\\ & p(w^T{x}_i+ϵ\ge 0)=p(ϵ\ge -w^T{x}_i)& \text{(9.104)}\\ & 1-\Phi (-w^T{x}_i)=\Phi (w^T{x}_i)& \text{(9.105)}\end{array}$$

上面用到了高斯分布的对称性(symmetry).这个潜在变量解释提供了你和模型的另外一种方法,具体参考本书11.4.6.

如果我们队$\delta$使用一个冈贝尔分布(Gumbel distribution),得到的就是对${ϵ}_i$的逻辑分布(logistic distribution),模型也降低到了(reduces to)逻辑回归(logistic regression).更多细节参考本书24.5.1.

9.4.3 有序概率回归(Ordinal probit regression)*

概率回归的潜在变量解释的一个优势就是很容易扩展到响应变量有序的情况下,也就是取C个离散值,以某种方式排序,比如从小到大等.这样就成了有序回归(ordinal regression).基本思想如下所述.引入$C+1$个阈值${\gamma }_j$以及集合:

$y_i=j\text{ if }{\gamma }_{j-1}\le z_i\le {\gamma }_j$(9.106)

其中的${\gamma }_0\le ...\le {\gamma }_C$.为了好辨认,就设${\gamma }_0=-\infty ,{\gamma }_1=0,{\gamma }_C=\infty$.例如,如果$C=2$,这就降低到了标准二值化概率模型(standard binary probit model),其中$z_i<0$导致$y_i=0$,而$z_i\ge 0$对应就是$y_i=1$.如果C=3,就把这个实线(real line)分成了三个区间:$(-\infty ,0],(0,{\gamma }_2],({\gamma }_2,\infty )$.这时候通过调节参数${\gamma }_2$就可以确保每个区间有正确相关规模的概率质量,也能符合每个类标签的经验频率(empirical frequencies).

对这个模型找最大似然估计(MLE)就可能比二值化概率回归要更麻烦点了,因为需要优化$w$和$\gamma$,而后面的那个必须服从一个有序约束(ordering constraint).(Kawakatsu and Largey 2009)提出了一个基于期望最大化(EM)的方法.另外从这个模型也可以推导出一个简单的吉布斯采样算法(Gibbs sampling algorithm),参考(Hoff 2009, p216).

9.4.4 多项概率模型(Multinomial probit models)*

然后考虑响应变量可以去C个无排序分类值的情况,也就是 y i ∈ { 1 , . . . , C } y_i\in \{1,...,C\} yi​∈{1,...,C}.多项概率模型定义如下所示:

$$\begin{array}{rlr}{z}_{ic}& =w^T{x}_{ic}+{ϵ}_{ic}& \text{(9.107)}\\ ϵ& \sim N(0,R)& \text{(9.108)}\\ y_i& =\arg \underset{c}{\max }{z}_{ic}& \text{(9.109)}\end{array}$$

关于这个模型的更多细节以及和多项逻辑回归(multinomial logistic regression)之间的关系可以去参考(Dow and Endersby 2004; Scott 2009; Fruhwirth-Schnatter and Fruhwirth 2010).(通过定义$w=[w_1,...,w_C],x_{ic}=[0,...,0,x_i,0,...,0]$,就可以恢复到更熟悉的函数方程$z_{ic}=x_i^T{w}_C$.)由于只有相对效用(relative utilities)有影响,可以将R限制为一个相关矩阵(correlation matrix).如果不去设$y_i=\arg {\max }_C{z}_{ic}$而是使用$y_{ic}=I(z_{ic}>0)$,就得到了多元概率模型(multivariate probit model),也是对于二值化输出相关的C建模的一种方法,参考(Talhouk et al. 2011).

9.5 多任务学习(multi-task learning)

有时候要对多个相关的分类或者拟合模型进行你和.通常都可能会假设不同模型之间的输入输出映射是相似的,所以可以在同时拟合所有参数来得到更好的性能.在机器学习里面,这一般叫做多任务学习 (multi-task learning,Caruana 1998), 转换学习 (transfer learning,Raina et al.2005), 或者(learning to learn,Thrun and Pratt 1997).在统计学上,通常的解决方法是分层贝叶斯模型(hierarchical Bayesian models,Bakker and Heskes 2003),当然也有一些其他方法(Chai 2010),后文会讲到.

9.5.1 多任务学习的分层贝叶斯

设$y_{ij}$是第j群(group)当中的第i项(item)的响应变量,其中$i=1:N_j,j=1:J$.例如,j可以对学校进行索引,然后i就是检索该学校中的学生,然后$y_{ij}$就是这个学生的测试分数,如本书5.6.2所示.或者也可以用j来对人进行索引,然后i代表队是购买次数,这样$y_{ij}$就只带被购买的特定商品(这就叫做离散选择模型(discrete choice modeling, Train 2009).设$x_{ij}$是对应$y_{ij}$的特征向量.目标就是要对于所有的j拟合出模型$p(y_j|x_j)$.

虽然有的组可能有很多数据,通常都是长尾的(long tail),其中大多数组都只有很少的数据.因此不能很有把握地去分开拟合各个模型,可又不想对所有组使用同一个模型.所以就折中一下,可以对每个组拟合一个离散的模型,但设置在不同的组织间模型参数具有相似性.更具体来说就是设$\mathrm{E}[y_{ij}|x_{ij}]=g(x_i^T{\beta }_j)$,其中的g是广义线性模型(GLM)的连接函数(link function).另外设${\beta }_j\sim N({\beta }_{\ast },{\sigma }_j^{2}I),{\beta }_{\ast }\sim N(\mu ,{\sigma }_{\ast }^{2}I)$.在这个模型里面,小样本规模的组就从更多样本的组借用统计强度,因为${\beta }_j$是和潜在通用母变量(latent common parents)${\beta }_{\ast }$相关的(这一点相关内容参考本书5.5).${\sigma }_j^2$这一项控制了第j组对通用母变量(common parents)的依赖程度,而${\sigma }_{\ast }^2$这一项控制了全局先验的强度.

为了简单起见,假设$\mu =0$,这样${\sigma }_j^2$和${\sigma }_{\ast }^2$就都是一直的了(可以通过交叉验证来设置).全局对数概率函数(overall log probability)形式为:

$\log p(D|\beta )+\log p(\beta )={\sum }_j[\log p((D_j|{\beta }_j)-\frac{||{\beta }_j-{\beta }_{\ast }|{|}^2}{2{\sigma }_j^2})]-\frac{||{\beta }_{\ast }|{|}^2}{2{\sigma }_{\ast }^2}$(9.110)

可以使用标准梯度方法进行对$\beta =({\beta }_{1:j},{\beta }_{\ast })$的最大后验估计(MAP).或者也可以使用迭代优化的策略,在${\beta }_j$和${\beta }_{\ast }$之间进行优化;因为似然函数和先验都是凸函数,这就保证了会收敛到全局最优解.要记住一旦一个模型训练出来了,就可以不用理会${\beta }_{\ast }$,分别使用每个模型.

9.5.2 应用:个性化垃圾邮件过滤

多任务学习的一个有趣应用就是个性化垃圾邮件过滤(personalized email spam filtering).假如要针对每个用户来拟合一个分类器${\beta }_j$.由于大部分用户都不会来标记他们的邮件是不是垃圾邮件,这就很难分开来对他们各自的模型进行估计.所以要设${\beta }_j$有一个通用的先验${\beta }_{\ast }$,表示了一个通用用户的参数.

这时候就可以利用上面的模型来估计这行为,这需要一点小技巧(Daume 2007b; Attenberg et al. 2009; Weinberger et al. 2009):将每个特征$x_i$都只做两个分本,一个联接(concatenated)到用户id,另外的一份则不做联接.这样要学习的预测器(predictor)的形式就为:

$\mathrm{E}[y_i|x_i,u]=({\beta }_{\ast },w_1,...,w_J{)}^T[x_i,I(u=1)x_i,...,I(u=J)x_i]$(9.111)

其中的u是用户id.也就是说:

$\mathrm{E}[y_i|x_i,u=j]=({\beta }_{\ast }+w_j{)}^T{x}_i$(9.112)

因此${\beta }_{\ast }$就可以从每个人的邮件中来进行估计,而$w_j$则只从第j个用户的邮件中来估计.

这和上面的分层贝叶斯模型具有对应关系(correspondence),定义$w_j={\beta }_j-{\beta }_{\ast }$.然后原始模型的对数概率函数就可以重新写成下面的形式:

${\sum }_j[\log p(D_j|{\beta }_{\ast }+w_j)-frac||w_j|{|}^{2}2{\sigma }_j^2]-\frac{||{\beta }_{\ast }|{|}^2}{2{\sigma }_{\ast }^2}$(9.113)

如果我们假设${\sigma }_j^2={\sigma }_{\ast }^2$,效果就和使用增强特征技巧(augmented feature trick)一样了,对$w_i$和${\beta }_{\ast }$都有一样的规范化强度(regularizer strength).不过如果让这两个不相等通常能得到更好的性能(Finkel and Manning 2009).

9.5.3 应用:域自适应(Domain adaptation)

域自适应问题是要训练从不同分布取样来的数据,比如邮件和新闻报道的文本.这个问题很明显是一个多任务学习的特例,其中各个任务都相同.

(Finkel and Manning 2009)使用了上面的分层贝叶斯模型来使用域自适应来实现两个自然语言处理任务(NLP tasks),命名实体识别(namely named entity recognition)和解译(parsing).他们的研究成果比区分每个数据集来拟合分散模型有更大的进步,而对于将所有数据放到一个简单模型来拟合相比就提升较小了.

9.5.4 其他类别的先验

在多任务学习里面,通常都假设先验是高斯分布的.不过有时候也可能选择别的先验更适合.比如对于联合分析(conjoint analysis)的任务,这个任务需要想办法弄清用户最喜欢产品的哪个特征.这可以使用跟前文同样的分层贝叶斯模型设置,不过要对${\beta }_j$使用更加稀疏的先验(sparsity-promoting prior),而不是高斯先验.这就叫做多任务特征选择(multi-task feature selection).更多可用方法等等参考(Lenk et al. 1996; Argyriou et al. 2008).

总去假设所有任务都一样想死并不总是很合理的.如果把不同质量(qualitatively different)的任务的参数汇集到一起,计算性能回避不汇聚要更差,因为咱们先验中的归纳偏见(inductive bias)是错误的.实际上已经有研究发现有点时候多任务学习要比分开解决每个任务效果更差,这也叫做负转换(negative transfer).

这类问题的另外一个方法就是使用更灵活的先验,比如混合高斯模型.这类灵活先验可以提供健壮性,应对先验误判(prior mis-specification).更多细节参考(Xue et al.2007; Jacob et al. 2008).当然也可以把稀疏提升先验(sparsity-promoting priors)混合起来使用(Ji et al. 2009).当然也有很多其他变体方法.

9.6 广义线性混合模型(Generalized linear mixed models)*

假如将多任务学习场景扩展来允许响应变量包含分组水平(group level)$x_j$,和项目水平(item level)$x_{ij}$.然后也允许参数${\beta }_j$在组间改变,或者参数$\alpha$在组间固定.这样就得到了下面的模型:

$\mathrm{E}[y_{ij}|x_{ij},x_j]=g(\varphi (x_{ij}{)}^T{\beta }_j+{\varphi }_2(x_j{)}^T{\beta }_j^{\prime }+{\varphi }_3(x_{ij}{)}^T\alpha +{\varphi }_t(x_j{)}^T{\alpha }^{\prime })$(9.114)
其中的${\varphi }_k$是基函数(basis functions).这个函数如图9.2(a)所示.(这种图在本书第十章会解释.)参数${\beta }_j$的个数随着组个数增长而增长,而参数$\alpha$则是固定的.

频率论里面将${\beta }_j$这一项叫做随机效应(random effects),因为它们在各组中随机变化;称$\alpha$为固定效应(fixed effect),看做是一个固定但未知的敞亮.如果一个模型同时有固定效应和随机效应,就称之为混合模型(mixed model).如果$p(y|x)$是一个广义线性模型(GLM),整个模型就叫做广义线性混合模型(generalized linear mixed effects model,缩写为GLMM).这种模型在统计学里面用的很普遍.

9.6.1 样例:针对医疗数据的半参数化广义线性混合模型(semi-parametric GLMMs for medical data)

下面的例子来自(Wand 2009).假如$y_{ij}$是第j个病人在第i次测量的脊椎骨矿物密度(spinal bone mineral density,缩写为SBMD).设$x_{ij}$为这个人的年龄,设$x_j$为此人的种族(ethnicity),可以是白种人/亚裔/黑种人/拉丁裔.主要目标就是要判断四个种族的平均脊椎骨矿物密度是否有显著区别.数据如图9.2(b)中浅灰色线所示.其中可见脊椎骨矿物密度(SBMD)和年龄有一个非线性关系,所以可以使用一个半参数化模型(semi-parametric model)综合了线性回归和非参数化回归(Ruppert et al. 2003).另外还可以发现每个组内不同个体之间也有变化,所以要用一个混合效应模型.具体来说就是使用${\varphi }_1(x_{ij})=1$来计量每个人的随机效应;${\varphi }_2(x_{ij})=0$,因为没有其他的按照人来变化的量;${\varphi }_3(x_{ij})=[b_k(x_{ij})]$,其中的$b_k$是第k个样条基函数(spline basis functions)(参考本书15.4.6.2),这是用来及计数年龄的非线性效应(nonlinear effect);另外还有个计算不同种族效应的${\varphi }_4(x_j)=[I(x_j=w),I(x_j=a),I(x_j=b),I(x_j=h)]$.然后使用一个线性连接函数(linear link function).整个模型就是:

$$\begin{array}{rlr}\mathrm{E}[y_{ij}|x_{ij},x_j]& ={\beta }_j+{\alpha }^{T}b(x_{ij})+{ϵ}_{ij}& \text{(9.115)}\\ & +{\alpha }_w^{\prime }I(x_j=w)+{\alpha }_a^{\prime }I(x_j=a)+{\alpha }_b^{\prime }I(x_j=b)+{\alpha }_h^{\prime }I(x_j=h)& \text{(9.116)}\end{array}$$

其中${ϵ}_{ij}\sim N(0,{\sigma }_y^2)$. $\alpha$包含了模型中和年龄相关的非参数部分,${\alpha }^{\prime }$包含了和种族相关的参数部分,而${\beta }_j$则包含了随着每个人j变化的随机偏移量.将这些回归系数(regression coefficients)都赋予各自的高斯先验.然后可以进行后验推导来计算$p(\alpha ,{\alpha }^{\prime },\beta ,{\sigma }^2|D)$(计算细节参考本书9.6.2).你喝了模型之后,可以对每个组计算预测.结果如图9.2(b)所示.还可以进行显著性检验(significance testing),以某个种群作为基准值(比如白人),来对每个种群g计算$p({\alpha }_g-{\alpha }_w|D)$,这就跟本书5.2.3里面讲的一样了.

此处参考原书图9.2

9.6.2 计算问题

广义线性混合模型(GLMMs)的主要问题是不太好拟合,这有两个原因.首先是$p(y_{ij}|\theta )$可能和先验$p(\theta )$未必共轭,其中$\theta =(\alpha ,\beta )$.另外一个原因是这个模型里面有两个未知层次,回归系数$\theta$以及先验$\eta =(\mu ,\sigma )$的均值和方差.

一个方法是采用全贝叶斯方法(fully Bayesian inference methods),比如变分贝叶斯(variational Bayes,Hall et al.2011),或者马尔科夫链蒙特卡罗方法(MCMC,Gelman and Hill 2007).变分贝叶斯(VB)在本书21.5会讲到,而马尔科夫链蒙特卡罗方法(MCMC)在本书24.1.

另一种方法就是使用经验贝叶斯(empirical Bayes),在本书5.6大概讲过.在广义线性混合模型的语境下,可以使用期望最大化算法(EM algorithm,本书11.4).
其中的E步骤计算$p(\theta |\eta ,D$,而M这一步骤优化$\eta$.如果是线性回归的情况下,E步骤就可以确切进行,但一般都不用确切计算出来,使用个近似值就行了.传统方法是使用数值正交(numerical quadrature)或者蒙特卡罗方法(Breslow and Clayton 1993).更快的方法是使用变分期望最大化(variational EM),参考Braun and McAuliffe 2010提供了对多水平离散选择模型问题使用变分期望最大化的应用案例.

在频率论统计学中,有一个拟合广义线性混合模型的流行方法叫做广义估计方程(generalized estimating equations,缩写为GEE,Hardin and Hilbe 2003).不过不推荐这个方法,因为在统计学上这个方法效率不如似然函数方法(参考本书6.4.3).另外这个方法也只能对人口参数$\alpha$进行估计,而不能对随机效应${\beta }_j$进行估计,而后者可能是更需要的.

$p(y_{ij}|\theta )$ $p(\theta )$ $\theta =(\alpha ,\beta )$ $\eta =(\mu ,\sigma )$

$p(\theta |\eta ,D)$

9.7 学习排序(Learning to rank)*

在本节降低是学习排序(learning to rank,缩写为LETOR)问题.也就是要学习一个函数,可以将一系列项目进行排序(后面会详细说具体内容).最常见的用途是信息检索(information retrieval).假如有一个检索查询(query)q,以及一系列文档$d^1,...,d^m$,这些文档和q可能相关(比如所有文档都包含字符串q).然后我们想要对文档按照和q的相关性来降序排列,展示出其中前面k个给用户.类似问题也出现在其他领域,比如协作过滤(collaborative filtering)(在一个游戏里面对玩家进行排名或者类似的问题,具体参考本书22.5.5).

接下来简单总结一些解决这个问题的方法,这部分内容参考了(Liu 2009).这部分内容其实并不是基于广义线性模型(GLM)的,但这本书其他地方也不适合放这些内容,就放这里了,就是这么任性.

衡量文档d和查询q之间相关性的标准方法是使用一个概率语言模型(probabilistic language model),基于词汇袋模型(bag of words model).也就是定义$sim(q,d)\stackrel{△}{=}p(q|d)={\prod }_{i=1}^{n}p(q_i|d)$,其中的$q_i$是第i个查询项目或者词汇,而$p(q_i|d)$是从文档d里面估计的一个多重伯努利分布(multinoulli distribution).实际使用中需要将这个估计分布进行光滑处理,比如使用狄利克雷先验(Dirichlet prior),来表示每个词汇的全局频率(overall frequency).这可以从系统中的全部文档来估计.具体来说就是使用下面这个表达式:

$p(t|d)=(1-\lambda )\frac{TF(td)}{LEN(d)}+\lambda p(t|background)$(9.117)

其中$TF(t|d)$是文档d中第t项的频率,而$LEN(d)$是文档d中词汇长度,而$0<\lambda <1$是光滑参数(smoothing parameter,更多细节参考 Zhai and Lafferty 2004).

不过也可能还有很多其他的信号也能用来衡量相关性.比如网络文本的页面评级(PageRank)就是对其权威性(authoritativeness)的衡量,是从网页链接结构来推导的(具体参考本书17.2.4).可以计算在一个文档中某个查询项目出现的次数和位置.接下来就讲一下如何将这些信号集中起来学习使用.

9.7.1 单点法(pointwise approach)

假如我们要收集代表一系列文档对每个查询项的相关度的训练数据.具体来说就是对每个查询q,设找到了m个可能的相关文档 $d_j,\text{for}j=1:m$.对每个查询文档对,定义一个特征向量$x(q,d)$.比如可能包含了查询项和文档的相似度排序,以及文档的页面评级分数(page rank score).然后假设有一系列的标签$y_i$代表的是文档$d_j$和查询项q之间的相关程度.这样的类标签可以十二指华东(比如是相关或不想管),或者也可以使用离散的相关程度来描述(比如很相关,有点相关,不相关).这样的类标签可以从查询项日志(query logs)获得,通过限制一个文档对于一个给定的查询项被点击的次数.

如果使用二值化相关标签,就可以使用标准二值化分类来估计$p(y=1|x(q,d))$,来解决这个问题你.如果使用排序相关标签,就可以使用有序回归(ordinal regression)$p(y=r|x(q,d))$来预测排序.这两种情况下都可以通过排序矩阵(scoring metric)来整理文档.这就叫做学习排序(LETOR)的单点法(pointwise approach),用的很广泛,因为特别简单.不过这个方法没有将文档在列表中的各自位置考虑进去.因此对列表中最末项目和最首位项目有一样的错误惩罚,这通常就不符合预期了.另外每次对相关性的判断也都非常短视(myopically).

9.7.2 成对法(pairwise approach)

Carterette et al. 2008 的研究证明人类要更擅长判断两个对象之间的相对相性,而不是绝对相关性.结果就导致数据可能告诉我们$d_j$比$d_k$和给定的查询更相关,或者反过来.可以对这种数据使用二值化分类器来进行建模,形式为$p(y_{jk}|x(q,d_j),x(q,d_k))$ ,当$rel(d_j,q)>rel(d_k,q)$时设$y_{jk}=1$,反之则设置$y_{jk}=0$.

对这个函数建模的一种方法如下所示:

$p(y_{jk}=1|x_j,x_k)=sigm(f(x_j)-f(x_k))$(9.118)

其中的$f(x)$是排序函数,一般都设置为线性函数$f(x)=w^{T}x$.这是一类特殊的神经网络,叫做排序网络(RankNet,Burges et al. 2005,关于神经网络的讨论参考本书16.5).可以通过最大化对数似然函数来找到w的最大似然估计(MLE),或者也可以等价地对交叉熵损失函数(cross entropy loss)最小化,也就是:

$$\begin{array}{rlr}L& =su{m}_{i=1}^N\sum _{j=1}^{m_i}\sum _{k=j+1}^{m_i}{L}_{ijk}& \text{(9.119)}\\ -L_{ijk}& =I(y_{ijk}=1)\log p(y_{ijk}=1|x_{ij},x_{ik},w)& \\ & +I(y_{ijk}=0)\log p(y_{ijk}=0|x_{ij},x_{ik},w)& \text{(9.120)}\end{array}$$

这就可以使用梯度下降法来优化了.微软的Bing搜索引擎就使用了排序网络的一个变种.

9.7.3 列表法(listwise approach)

承兑发的问题就是关于相关性的决策判断只是基于一对项目或者一堆文档,而不是考虑完整的语境(context).接下来要讲的方法就要同时查看整个项目列表.

在列表上可以制定一个索引排序来定义一个全排列,$\pi$.要对关于$\pi$的不确定性建模,就可以使用一个Plackett-Luce分布,这个分布的命名就是基于两个独立推导该公式的人(Plackett 1975)和(Luce 1959).这个分布形式如下所示:

$p(\pi |s)={\prod }_{j=1}^m\frac{s_j}{{\sum }_{u=j}^m{s}_u}$(9.121)

其中的$s_j=s({\pi }^{-1}(j))$是对排在第j个位置的文档的排序.

要理解等式9.121,可以举个简单例子.设$\pi =(A,B,C)$.然后就得到了$p(\pi )$是A排第一的概率,乘以A排第一的条件下B排第二的概率,再乘以AB分别排第一第二之后C排第三未知的概率,也就是说:

$p(\pi |s)=\frac{s_A}{s_A+s_B+s_C}\times \frac{s_B}{s_B+s_C}\times \frac{s_C}{s_C}$(9.122)

要整合特征,可以定义$s(d)=f(x(q,d))$,通常都设f为线性函数$f(x)=w^{T}x$.这就叫做列表网络模型(ListNet model,Cao et al. 2007).要训练这个模型,设$y_i$是文档对于查询项目i的相关性分数.然后就要最小化交叉熵项(cross entropy term):

$-{\sum }_i{\sum }_{p}ip(\pi |y_i)\log p(\pi |s_i)$(9.123)

当然了,这也很困难,因为第i项目需要累加总共超过$m_i!$次排列才行.要让这个好处理,可以考虑只在前面的k个位置上进行排列:

$p({\pi }_{1:k}|s_{1:m})={\prod }_{j=1}^k\frac{s_j}{{\sum }_{u=1}^m{s}_u}$(9.124)

这样就只需要$m!/(m-k)!$次这样的排列了.如果设置$k=1$,那就可以在$O(m)$时间内计算每个交叉熵以及其导数了.

在列表中只有一个文档被认为相关的特例下,比如设$y_i=c$,就可以使用多项逻辑回归(multinomial logistic regressi)了:

$p(y_i=c|x)=\frac{\exp (s_C)}{{\sum }_{(}^m{c}^{\prime }=1)\exp (s_{c^{\prime }})}$(9.125)

这个方法至少能达到排序方法的性能水平,至少在协作过滤(collaborative filtering)的情况下如此(Yang et al. 2011).

9.7.4 排序的损失函数

对一个排序系统的性能衡量,有好几种方法,大概如下所述.

• 平均排序倒数(Mean reciprocal rank,缩写为MRR).对于一个查询项q,设其第一个相关文档的排序位置记作$r(q)$.然后定义一个平均排序倒数为$1/r(q)$.这是很简单的性能衡量.
• 均值平均准确率(Mean average precision,缩写为MAP,注意要和最大后验分布 MAP相区分).在二值化相关标签的情况下,可以定义某个排序在k位置上的精度如下:

$P@k(\pi )\stackrel{△}{=}\frac{\text{num. relevant documents in the top k positions of}\pi }{k}$(9.126)

然后可以定义平均精度(average precision)如下:

$AP(\pi )\stackrel{△}{=}\frac{{\sum }_{k}P@k(\pi )\times I_k}{\text{num relevant documents}}$(9.127)

其中当且仅当文档k为相关的时候$I_k$才等于1.例如,如果有相关标签$y=(1,0,1,0,1)$,然后AP就是 $\frac{1}{3}(\frac{1}{1}+\frac{2}{3}+\frac{3}{5})\approx 0.76$最终就定义了均值平均精度(mean average precision)为对所有查询上的AP求平均值.

• 归一化折扣累积增益(Normalized discounted cumulative gain,缩写为NDCG).假如香瓜鸟枪有多种层次.就可以定义对前面以一定次序排列的k个项目的折扣累积增益(discounted cumulative gain)如下所示:

$DCG@k(r)=r_1+{\sum }_{i=2}^k\frac{r_i}{{\log }_{2}i}$(9.128)

其中的$r_i$是第i项的相关性,而${\log }_2$项目是用来稍后在列表中扣除项目的.表9.3展示了一个简单的数值样本.另一重定义就是强调了检索到相关文档(retrieving relevant documents),使用的是:

$DCG@k(r)={\sum }_{i=1}^k\frac{2^{r_i}-1}{{\log }_2(1+i)}$(9.129)

折扣累积增益(DCG)的一个问题就是只要返回列表的长度变化,这个增益的数量级就会有变化.因此通常都要用理想折扣累积增益(ideal DCG)来将其归一化,理想折扣累积增益(ideal DCG)是指通过使用最优排序来得到的DCG,即$IDCG@k(r)=\arg {\max }_{\pi }DCG@k(r)$.最终就定义出来了归一化折扣累积增益(Normalized discounted cumulative gain,缩写为NDCG),定义为$DCG/IDCG$.表9.4给出了一个简单数值样本.NDCG方法可以对查询项目进行平均然后来衡量性能.

• 排序相关性(Rank correlation).可以在排序列表$\pi$和相关性判断${\pi }^{\ast }$之间衡量相关性,使用的方法就很多了.比如可以使用加权肯德尔$\tau$统计(weighted Kendall’s τ statistics),这个统计定义形式为两个列表间不连续的加权值对:

$\tau (\pi ,{\pi }^{\ast })=\frac{{\sum }_{u<v}{w}_{uv}[1+sgn({\pi }_u-{\pi }_v)sgn({\pi }_u^{\ast }-{\pi }_v^{\ast })]}{2{\sum }_{u<v}{w}_{uv}}$(9.130)

其他方法也都很常用.

这些损失函数可以有不同用法.在贝叶斯方法中,首先使用后验推断来拟合模型;这就一来似然函数和先验,但不用管损失函数.然后选在测试的时候选择行为来最小化期望未来损失(expected future loss).一种方法就是从后严重对参数取样,即${\theta }^s\sim p(\theta |D)$,然后评估,比如对在k精确度不同的阈值,在${\theta }^s$上取平均值.这种方法的具体样例参考(Zhang et al. 2010).

在频率论方法中,在训练集上就要试图最小化经验损失函数.问题就是这些损失函数并不是模型参数的可微函数(differentiable functions).要么就要用非梯度的优化方法,要么就要使用代理损失函数来替代进行最小化.交叉熵损失函数(比如负对数似然函数)是一个广泛使用的代理损失函数.

另外一种损失函数,也叫作加权估计-排序成对损失函数(weighted approximate-rank pairwise loss,缩写为WARP loss),由(Usunier et al. 2009)提出,然后由(Weston et al. 2010)对其进行了扩展,提供了对在k精度上的损失函数的更好的估计.这个损失函数的定义如下所示:

$$\begin{array}{rlr}WARP(f(x,:),y)& \stackrel{△}{=}L(rank(f(x,:),y)& \text{(9.131)}\\ rank(f(x,:),y)& =\sum _{y^{\prime }\ne y}I(fx,y^{\prime })\ge f(x,y))& \text{(9.132)}\\ L(k)& \stackrel{△}{=}\sum _{j=1}^k{\alpha }_j,\text{with}{\alpha }_1\ge {\alpha }_{\ge }...\ge 0& \text{(9.133)}\end{array}$$

上式中的$f(x,:)=[f(x,1),...,f(x,|y|)]$是对每个可能输出标签的评分向量,或者在迭代重加权(IR)项目中,就是对每个对应输入查询x的每个可能文章的评分向量.表达式$rank(f(x,:),y)$衡量由该评分函数分配真实标签y的评分.最后的L是讲整数值的评分转化成实数值的惩罚项(real-valued penalty).使用${\alpha }_1=1,{\alpha }_{j>1}=0$可以优化前部位置分类标签正确的比例.设置${\alpha }_{1:k}$为非零值可以优化排序列表中的前k个项目,因为使用了均值平均准确率(MAP)或k精确度来衡量,所以性能不错.即便这样,加权估计-排序成对损失函数(weighted approximate-rank pairwise loss,缩写为WARP loss)也还是很难去优化的,但是可以使用蒙特卡罗取样方法来近似,然后再用梯度下降法去优化,这部分参考(Weston et al. 2010).

练习略