简述
在 Liner Regression 中 y \boldsymbol{y} y 是一个连续值,那么当我们解决一个二分类问题时 y \boldsymbol{y} y 则是一个离散值,只有两个取值( 0 0 0 或 1 1 1 ),这时可以通过广义线性模型来解决。通过广义线性模型,我们找到一个单调可微函数将分类任务的标记 y \boldsymbol{y} y 与线性模型的预测值联系起来。最后我们找到 Logistic Function 来作为线性模型。
广义线性模型(Generalized Linear Models )
在笔记中,我先把广义线性模型梳理一遍,这有助于我更好得学习机器学习。
当然如果仅仅是为了学习 Logistic Function 那么可以先看看后面的内容,之后再看这部分内容。
在回归学习中,我们的函数都类似于f ( y ∣ x ) ; θ ∼ N ( μ , σ 2 ) f(\boldsymbol{y}|\boldsymbol{x}); \boldsymbol{\theta} \sim \mathcal{N}(\mu,\sigma^2) f ( y ∣ x ) ; θ ∼ N ( μ , σ 2 ) 或者在之后讲的二分类函数 f ( y ∣ x ) ; θ ∼ B e r n o u l l i ( ϕ ) f(\boldsymbol{y}|\boldsymbol{x}); \boldsymbol{\theta} \sim \rm{Bernoulli}(\phi) f ( y ∣ x ) ; θ ∼ B e r n o u l l i ( ϕ ) ,这里的 μ \mu μ 和 ϕ \phi ϕ 都分别是 x x x 和 θ \theta θ 的某种函数( σ \sigma σ 与分布无关 )。其实,有一种更广泛的模型,这两种模型都是它的特例,这种更广泛的模型叫做广义线性模型。
在广义线性模型中(GLM), 对于每个独立参数的 y \boldsymbol{y} y ,假设通过一个指数族产生。这就是说,对于均值 μ \mu μ , 和独立变量 x \boldsymbol{x} x ,有:
E ( y ) = μ = g − 1 ( θ T x ) E(\boldsymbol{y})=\boldsymbol{\mu}=g^{-1}(\boldsymbol{\theta}^{\mathit{T}} \boldsymbol{x})
E ( y ) = μ = g − 1 ( θ T x )
E ( y ) E(\boldsymbol{y}) E ( y ) 是 y \boldsymbol{y} y 的期望;θ T x \boldsymbol{\theta}^{\mathit{T}} \boldsymbol{x} θ T x 是一个线性估计; g g g 是链接函数。
关于广义线性模型更多的知识请前往 Wikipedia 。
指数族(The Exponential Family )
在学习广义线性模型之前,我们要先定义一下指数族分布(exponential family distributions)。如果一个分布能用下面的方式来写出来,我们就说这类分布属于指数族:
p ( y ; η ) = b ( y ) e x p ( η T T ( y ) − a ( η ) ) p(y;\eta)=b(y)\rm{exp}(\eta^TT(y)-a(\eta))
p ( y ; η ) = b ( y ) e x p ( η T T ( y ) − a ( η ) )
上面的式子中,
η \eta η :该分布的自然参数(natural parameter,也叫典范参数 canonical parameter);
T ( y ) T(y) T ( y ) :充分统计量(sufficient statistic),我们目前用的这些分布中通常 T ( y ) = y T(y) = y T ( y ) = y ;
a ( η ) a(\eta) a ( η ) :一个对数分割函数(log partition function) ;
$e^{−a(\eta)} $ :这个量本质上扮演了归一化常数(normalization constant)的角色,也就是确保分布的 p ( y ; η ) p(y;\eta) p ( y ; η ) 的总和等于1。
对 给定的 T T T , a a a 和 b b b 就定义了一个以 η \eta η 为参数的分布族(family,或者叫集 set);通过改变 η \eta η ,我们就能得到这个分布族中的不同分布。
现在咱们看到的伯努利(Bernoulli)分布和高斯(Gaussian)分布就都属于指数分布族。伯努利分布的均值是 ϕ \phi ϕ ,也写作 B e r n o u l l i ( ϕ ) \rm{Bernoulli}(\phi) B e r n o u l l i ( ϕ ) ,确定的分布是 y ∈ { 0 , 1 } y \in\{0,1\} y ∈ { 0 , 1 } ,因此有 p ( y = 1 ; ϕ ) = ϕ ; p ( y = 0 ; ϕ ) = 1 − ϕ p(y=1;\phi)=\phi;p(y=0;\phi)=1-\phi p ( y = 1 ; ϕ ) = ϕ ; p ( y = 0 ; ϕ ) = 1 − ϕ 。这时候只要修改 ϕ \phi ϕ ,就能得到一系列不同均值的伯努利分布了。现在我们展示的通过修改 ϕ \phi ϕ ,而得到的这种伯努利分布,就属于指数分布族;也就是说,只要给定一组 T T T , a a a 和 b b b ,就可以用上面的等式来确定一组特定的伯努利分布了。
伯努利分布通过广义线性模型可以这样写:
p ( y ; ϕ ) = ϕ y ( 1 − ϕ ) 1 − y = e x p ( y log ϕ + ( 1 − y ) log ( 1 − ϕ ) ) = e x p ( ( log ( ϕ 1 − ϕ ) ) y + log ( 1 − ϕ ) ) \begin{aligned}
p(y;\phi) &= \phi^y(1-\phi)^{1-y}\\
&=\rm{exp}(\it{y}\rm{}\log{\phi}+(1-y)\log{(1-\phi)})\\
&=\rm{exp}\left(\left(\log{\left(\frac{\phi}{1-\phi}\right)}\right)\it{y}\rm{}+\log{(1-\phi)}\right)
\end{aligned}
p ( y ; ϕ ) = ϕ y ( 1 − ϕ ) 1 − y = e x p ( y log ϕ + ( 1 − y ) log ( 1 − ϕ ) ) = e x p ( ( log ( 1 − ϕ ϕ ) ) y + log ( 1 − ϕ ) )
因此,给出了自然参数(natural parameter)η = log ( ϕ 1 − ϕ ) \eta=\log{\left(\frac{\phi}{1-\phi}\right)} η = log ( 1 − ϕ ϕ ) 。 有趣的是,如果我们翻转这个定义,通过 ϕ \phi ϕ 表示 η \eta η 就会得到 ϕ = 1 1 + e − η \phi=\frac{1}{1+e^{-\eta}} ϕ = 1 + e − η 1 。这正好就是我们之后在 Logistic Function 中会见到的 S 型函数(sigmoid function)! 在我们把逻辑回归作为一种广义线性模型(GLM)的时候还会遇到伯努利分布以如下情况表示。
T ( y ) = y a ( η ) = − log ( 1 − ϕ ) = l o g ( 1 + e η ) b ( y ) = 1 \begin{aligned}
T(y) &= y\\
a(\eta)&=-\log{(1-\phi)}\\
&=log(1+e^\eta)\\
b(y)&=1
\end{aligned}
T ( y ) a ( η ) b ( y ) = y = − log ( 1 − ϕ ) = l o g ( 1 + e η ) = 1
接下来就看看高斯分布。在推导线性回归的时候, σ 2 \sigma^2 σ 2 的值对我们最终选择的 $\theta $ 和 h θ ( x ) h_\theta(x) h θ ( x ) 都没有影响。所以我们可以给 σ 2 \sigma^2 σ 2 取一个任意值。为了简化推导过程,就令 σ 2 = 1 \sigma^2 = 1 σ 2 = 1 。然后就有了下面的等式:
p ( y ; μ ) = 1 2 π e x p ( − 1 2 ( y − μ ) ) = 1 2 π e x p { − 1 2 y 2 } ⋅ e x p { μ y − 1 2 μ 2 } \begin{aligned}
p(y;\mu)&=\frac{1}{\sqrt{2\pi}}\rm{exp}\left(-\frac{1}{2}(\mathcal{y}\rm{}-\mu)\right)\\
&=\frac{1}{\sqrt{2\pi}}\rm{exp} \left \{ {-\frac{1}{2}\mathcal{y}^2} \right \} \cdot \rm{exp}\left \{{\mu\mathcal{y}-\frac{1}{2}\mu^2} \right \}
\end{aligned}
p ( y ; μ ) = 2 π 1 e x p ( − 2 1 ( y − μ ) ) = 2 π 1 e x p { − 2 1 y 2 } ⋅ e x p { μ y − 2 1 μ 2 }
注:如果我们把 σ 2 \sigma^2 σ 2 作为一个变量,高斯分布就也可以表达成指数分布的形式,其中 η ∈ R 2 \eta \in \mathbb{R}^2 η ∈ R 2 就是一个同时依赖 μ \mu μ 和 σ \sigma σ 的二维向量。然而,对于广义线性模型GLMs方面的用途,参数 σ 2 \sigma^2 σ 2 就也可以看成是对指数分布族的更泛化的定义:p ( y ; η , τ ) = b ( a , τ ) e x p ( ( η T T ( y ) − a ( η ) ) / c ( τ ) ) p(y;\eta,\tau)=b(a,\tau)\rm{exp}((\eta^T\it{T}(\mathcal{y})-\mathcal{a}(\eta))/\mathcal{c}(\tau)) p ( y ; η , τ ) = b ( a , τ ) e x p ( ( η T T ( y ) − a ( η ) ) / c ( τ ) ) 。这里面的 τ \tau τ 叫做分散度参数(dispersion parameter),对于高斯分布,c ( τ ) = σ 2 c(\tau)=\sigma^2 c ( τ ) = σ 2 ;不过上文中已经进行了简化,所以就不对各种需要考虑的情况进行更为泛化的定义了。
Logistic Regression
基本过程
从最基本开始,我们不使用广义线性模型,对于二分类问题,输出标记 $\boldsymbol{y} \in {0,1} $ ,于是我们使用线性回归最基本的模型 z = θ T x z= \boldsymbol{\theta}^{\mathit{T}} \boldsymbol{x} z = θ T x 来预测 y \boldsymbol{y} y ,即我们将 z z z 转化为 0 / 1 0/1 0 / 1 值。最理想的是 “单位越阶函数”(unit-step function):
y = { 0 , z < 0 ; 0.5 , z = 0 ; 1 , z > 0 ; y=\left \{
\begin{aligned}
0,& &z<0; \\
0.5,& &z=0; \\
1, & &z>0;
\end{aligned}
\right .
y = ⎩ ⎪ ⎨ ⎪ ⎧ 0 , 0 . 5 , 1 , z < 0 ; z = 0 ; z > 0 ;
但是,很明显,此函数不是很完美,于是,我们找到了一个“替代函数”来近似这个“单位跃阶函数”,并希望它单调可微,对数几率函数(Logistic Function)便满足这样一个条件:
g ( z ) = 1 1 + e − z g(z)=\frac{1}{1+e^{-z}}
g ( z ) = 1 + e − z 1
由图你能直观得看到,当 z → + ∞ z\to+\infty z → + ∞ 的时候 g ( z ) g(z) g ( z ) 趋向于 1 1 1 ,而当 z → − ∞ z\to-\infty z → − ∞ 时 g ( z ) g(z) g ( z ) 趋向于 0 0 0 。
其实 g ( z ) g(z) g ( z ) 也是 h θ ( x ) h_\theta(x) h θ ( x ) ,且,像最开始一样,我们规定 X 0 = 1 X_0 =1 X 0 = 1 ,于是有:θ T x = θ 0 + ∑ j = 1 n θ j x j \boldsymbol{\theta}^Tx=\theta_0+\sum^{n}_{j=1}\theta_jx_j θ T x = θ 0 + ∑ j = 1 n θ j x j
现在我们看看 g ′ ( z ) g'(z) g ′ ( z ) 的特性:
g ′ ( z ) = d d z 1 1 + e − z = 1 1 + e − z ( e − z ) = 1 1 + e − z ( 1 − 1 1 + e − z ) = g ( z ) ( 1 − g ( z ) ) . \begin{aligned}
g'(z) &= \frac{d}{dz}\frac{1}{1+e^{-z}}\\
&=\frac{1}{1+e^{-z}}(e^{-z})\\
&=\frac{1}{1+e^{-z}}\left(1-\frac{1}{1+e^{-z}}\right)\\
&=g(z)(1-g(z)).
\end{aligned}
g ′ ( z ) = d z d 1 + e − z 1 = 1 + e − z 1 ( e − z ) = 1 + e − z 1 ( 1 − 1 + e − z 1 ) = g ( z ) ( 1 − g ( z ) ) .
接着我们通过对 h θ ( x ) h_\theta(x) h θ ( x ) 进行假设,得到:
p ( y ∣ x ; θ ) = ( h θ ( x ) ) y ( 1 − h θ ( x ) ) 1 − y p(y|x;\theta)=(h_\theta(x))^y(1-h_\theta(x))^{1-y}
p ( y ∣ x ; θ ) = ( h θ ( x ) ) y ( 1 − h θ ( x ) ) 1 − y
然后就能写出似然函数:
L ( θ ) = p ( y ∣ X ; θ ) L(\boldsymbol{\theta}) = p(\boldsymbol{y}|X;\boldsymbol{\theta})
L ( θ ) = p ( y ∣ X ; θ )
又与之前一样写出 L ( θ ) L(\theta) L ( θ ) 的对数函数 ℓ ( θ ) \ell(\boldsymbol{\theta}) ℓ ( θ ) 以方便计算。
于是有 cost function:
J ( θ ) = − 1 m ℓ ( θ ) J(\theta) = -\frac{1}{m}\ell(\boldsymbol{\theta})
J ( θ ) = − m 1 ℓ ( θ )
然后目标就是:
θ = arg min J ( θ ) \boldsymbol{\theta} = \mathop{\argmin}_{}{J(\theta)}
θ = a r g m i n J ( θ )
θ \boldsymbol{\theta} θ 的求法
我们从最开始得到的假设函数讲起。如何得到它呢?
假设函数
我们首先假设:
P ( y = 1 ∣ x ; θ ) = h θ ( x ) P ( y = 0 ∣ x ; θ ) = 1 − h θ ( x ) P(y=1|x;\theta)=h_\theta(x)\\
P(y=0|x;\theta)=1-h_\theta(x)
P ( y = 1 ∣ x ; θ ) = h θ ( x ) P ( y = 0 ∣ x ; θ ) = 1 − h θ ( x )
于是,更简单的写法就是:
p ( y ∣ x ; θ ) = ( h θ ( x ) ) y ( 1 − h θ ( x ) ) 1 − y p(y|x;\theta)=(h_\theta(x))^y(1-h_\theta(x))^{1-y}
p ( y ∣ x ; θ ) = ( h θ ( x ) ) y ( 1 − h θ ( x ) ) 1 − y
似然函数(likelihood function)
假设 m m m 个训练样本都是各自独立的,那么就可以按如下的方式来写带参数的似然函数:
L ( θ ) = p ( y ∣ X ; θ ) = ∏ i = 1 m p ( y ( i ) ∣ x ( i ) ; θ ) 将 不 同 的 样 本 的 概 率 相 乘 = ∏ i = 1 m ( h θ ( x ( i ) ) ) y ( i ) ( 1 − h θ ( x ( i ) ) ) 1 − y ( i ) \begin{aligned}
L(\boldsymbol{\theta}) &= p(\boldsymbol{y}|X;\boldsymbol{\theta})\\
&=\prod^m_{i=1}p(y^{(i)}|x^{(i)};\boldsymbol{\theta}) &将不同的样本的概率相乘\\
&=\prod^m_{i=1}(h_\theta(x^{(i)}))^{y^{(i)}}(1-h_\theta(x^{(i)}))^{1-y^{(i)}}
\end{aligned}
L ( θ ) = p ( y ∣ X ; θ ) = i = 1 ∏ m p ( y ( i ) ∣ x ( i ) ; θ ) = i = 1 ∏ m ( h θ ( x ( i ) ) ) y ( i ) ( 1 − h θ ( x ( i ) ) ) 1 − y ( i ) 将 不 同 的 样 本 的 概 率 相 乘
取对数更容易计算:
ℓ ( θ ) = log L ( θ ) = ∑ i = 1 m y i log h ( x ( i ) ) + ( 1 − y i ) log ( 1 − h ( x ( i ) ) ) \begin{aligned}
\ell(\boldsymbol{\theta})&=\log{L(\theta)}\\
&=\sum^{m}_{i=1}y^{i}\log{h(x^{(i)})}+ (1-y^{i})\log{(1-h(x^{(i)}))}
\end{aligned}
ℓ ( θ ) = log L ( θ ) = i = 1 ∑ m y i log h ( x ( i ) ) + ( 1 − y i ) log ( 1 − h ( x ( i ) ) )
Cost Function
极大似然函数中,为了求得最优的 θ \boldsymbol{\theta} θ ,就是让 ℓ ( θ ) \ell(\boldsymbol{\theta}) ℓ ( θ ) 尽可能得大,于是在 cost function 中,我们加入系数 − 1 m -\frac{1}{m} − m 1 ,于是得到了:
J ( θ ) = − 1 m ∑ i = 1 m y i log h ( x ( i ) ) + ( 1 − y i ) log ( 1 − h ( x ( i ) ) ) J(\theta) = -\frac{1}{m}\sum^{m}_{i=1}y^{i}\log{h(x^{(i)})}+ (1-y^{i})\log{(1-h(x^{(i)}))}
J ( θ ) = − m 1 i = 1 ∑ m y i log h ( x ( i ) ) + ( 1 − y i ) log ( 1 − h ( x ( i ) ) )
梯度下降法(gradient descent)
其实,可以直接对 ℓ ( θ ) \ell(\boldsymbol{\theta}) ℓ ( θ ) 做梯度上升求得 θ \boldsymbol{\theta} θ 。
按照向量的形式,我们对 θ \theta θ 的更新可以写成:
θ : = θ − α ⋅ ∇ θ J ( θ ) \boldsymbol{\theta}:=\boldsymbol{\theta}-\alpha \cdot \nabla_{\theta} J(\boldsymbol{\theta})
θ : = θ − α ⋅ ∇ θ J ( θ )
找到最优的第一步是对 J ( θ ) J(\boldsymbol{\theta}) J ( θ ) 求导,我们从一个样本开始:
∂ ∂ θ j J ( θ ) = − 1 m ℓ ( θ ) = − 1 m ( y 1 g ( θ T x ) − ( 1 − y ) 1 1 − g ( θ T x ) ) ∂ ∂ θ j g ( θ T x ) = − 1 m ( y 1 g ( θ T x ) − ( 1 − y ) 1 1 − g ( θ T x ) ) g ( θ T x ) ( 1 − g ( θ T x ) ) = − 1 m ( y ( 1 − g ( θ T x ) ) − ( 1 − y ) g ( θ T x ) ) x j = 1 m ( h θ ( x ) − y ) x j \begin{aligned}
\frac{\partial}{\partial \theta_j}J(\theta)&=-\frac{1}{m}\ell(\theta)\\
&=-\frac{1}{m}\left( {y\frac{1}{g(\theta^Tx)}-(1-y)\frac{1}{1-g(\theta^Tx)}}\right)\frac{\partial}{\partial\theta_j}g(\theta^Tx)\\
&=-\frac{1}{m}\left( {y\frac{1}{g(\theta^Tx)}-(1-y)\frac{1}{1-g(\theta^Tx)}}\right)g(\theta^Tx)(1-g(\theta^Tx))\\
&=-\frac{1}{m}\left(y(1-g(\theta^Tx))-(1-y)g(\theta^Tx)\right)x_j\\
&=\frac{1}{m}(h_\theta(x)-y)x_j
\end{aligned}
∂ θ j ∂ J ( θ ) = − m 1 ℓ ( θ ) = − m 1 ( y g ( θ T x ) 1 − ( 1 − y ) 1 − g ( θ T x ) 1 ) ∂ θ j ∂ g ( θ T x ) = − m 1 ( y g ( θ T x ) 1 − ( 1 − y ) 1 − g ( θ T x ) 1 ) g ( θ T x ) ( 1 − g ( θ T x ) ) = − m 1 ( y ( 1 − g ( θ T x ) ) − ( 1 − y ) g ( θ T x ) ) x j = m 1 ( h θ ( x ) − y ) x j
其中,用到了上面提到的 g ′ ( z ) g'(z) g ′ ( z ) 的特性,即 ∂ ∂ θ j g ( θ T x ) = g ( θ T x ) ( 1 − g ( θ T x ) ) \frac{\partial}{\partial\theta_j}g(\theta^Tx)=g(\theta^Tx)(1-g(\theta^Tx)) ∂ θ j ∂ g ( θ T x ) = g ( θ T x ) ( 1 − g ( θ T x ) ) ,然后梯度上升就简单写为:
θ : = θ − α ⋅ 1 m ( y − h θ ( x ) ) x j \boldsymbol{\theta}:=\boldsymbol{\theta}-\alpha \cdot \frac{1}{m} (y-h_\theta(x))x_j
θ : = θ − α ⋅ m 1 ( y − h θ ( x ) ) x j
然后,再扩展为一个训练集:
Repeat until convergence { θ j : = θ j + α ⋅ 1 m ∑ i = 1 m ( y ( i ) − h θ ( x ( i ) ) ) x j } \begin{aligned}
&\text{Repeat until convergence}\{\\
& \quad \quad \theta_j := \theta_j +\alpha \cdot \frac{1}{m} \sum^{m}_{i=1}(y^{(i)}-h_{\theta}(x^{(i)}))x_j\\
&\}
\end{aligned}
Repeat until convergence { θ j : = θ j + α ⋅ m 1 i = 1 ∑ m ( y ( i ) − h θ ( x ( i ) ) ) x j }
有趣的是,这个式子正好与线性回归看上去一样,但是这实际上并不相同,原因是,我们对于 h θ ( x ) h_\theta(x) h θ ( x ) 的定义不同。但为什么相似呢?深层次的原因在于 广义线性模型 。
L ( θ ) L(\theta) L ( θ ) 最大的其它算法
下面这个方法更好,但是数学难度较高,其基本方法是“求方程零点的牛顿法”。
具体讲讲。假如我们有一个从实数到实数映射的函数 f : R ↦ R \it{f}:\rm{R} \mapsto \rm{R} f : R ↦ R ,然后要找到一个 θ \theta θ ,来满足 f ( θ ) = 0 \it{f}\,\rm{ }(\theta)=0 f ( θ ) = 0 ,其中θ ∈ R \theta \in R θ ∈ R 是一个实数。牛顿法就是对 θ \theta θ 进行如下的更新:
θ : = θ − f ( θ ) f ′ ( θ ) \theta:=\theta-\frac{f(\theta)}{f'(\theta)}
θ : = θ − f ′ ( θ ) f ( θ )
对这个公式的简单解释是:通过一条逼近曲线的直线(切线)不断迭代来找到零点。
对于上述的图,在 A 的切线方程为:y A = f ( x A ) + f ′ ( x A ) ( x − x A ) y_A = f(x_A)+f'(x_A)(x-x_A) y A = f ( x A ) + f ′ ( x A ) ( x − x A ) 。
为了求下一个迭代点 B ,则有 y B = 0 = f ( x A ) + f ′ ( x A ) ( x B − x A ) y_B = 0 = f(x_A)+f'(x_A)(x_B-x_A) y B = 0 = f ( x A ) + f ′ ( x A ) ( x B − x A ) ,
即:
⟹ x B = x a − f ( x A ) f ′ ( x A ) ⟹ x n + 1 = x n − f ( x n ) f ′ ( x n ) ⟹ θ : = θ − f ( θ ) f ′ ( θ ) \begin{aligned}
&\Longrightarrow x_B=x_a-\frac{f(x_A)}{f'(x_A)}\\\\
&\Longrightarrow x_{n+1}=x_{n}-\frac{f(x_n)}{f'(x_n)}\\\\
&\Longrightarrow \theta:=\theta-\frac{f(\theta)}{f'(\theta)}
\end{aligned}
⟹ x B = x a − f ′ ( x A ) f ( x A ) ⟹ x n + 1 = x n − f ′ ( x n ) f ( x n ) ⟹ θ : = θ − f ′ ( θ ) f ( θ )
用这个办法,我们求解 ℓ ( θ ) \ell(\theta) ℓ ( θ ) 中 θ \theta θ 的最优取值,即算 ℓ ( θ ) = 0 \ell(\theta) = 0 ℓ ( θ ) = 0 的解,即可以通过以下迭代式求得:
θ : = θ − ℓ ′ ( θ ) ℓ ′ ′ ( θ ) \theta:=\theta-\frac{\ell'(\theta)}{\ell''(\theta)}
θ : = θ − ℓ ′ ′ ( θ ) ℓ ′ ( θ )
对于向量 θ \boldsymbol{\theta} θ 的求解,我们扩展牛顿法到多维的情况,叫做牛顿-拉普森法(Newton-Raphson method),如下:
θ : = θ − H − 1 ∇ θ ℓ ( θ ) \boldsymbol{\theta} := \boldsymbol{\theta}-H^{-1}\nabla_{\theta} \ell(\boldsymbol{\theta})
θ : = θ − H − 1 ∇ θ ℓ ( θ )
其中 ∇ θ ℓ ( θ ) \nabla_{\theta} \ell(\boldsymbol{\theta}) ∇ θ ℓ ( θ ) 是 ℓ ( θ ) \ell(\boldsymbol{\theta}) ℓ ( θ ) 对 θ \theta θ 的偏导。H H H 是一个 n × n n\times n n × n 矩阵(考虑 θ 0 \theta_0 θ 0 的话是 ( n + 1 ) × ( n + 1 ) (n+1)\times (n+1) ( n + 1 ) × ( n + 1 ) ),也可叫做 Hessian,具体定义为:
H i j = ∂ 2 ℓ ( θ ) ∂ θ i ∂ θ j H_{ij} = \frac{\partial^2\ell(\theta)}{\partial\theta_i\partial\theta_j}
H i j = ∂ θ i ∂ θ j ∂ 2 ℓ ( θ )
注意,当 n n n 小的时候牛顿法的速度明显更快,但是当 n n n 较大时,由于需要处理 Hession 矩阵,时间开销急剧增加。
构建广义线性模型
设想你要构建一个模型,来估计在给定的某个小时内来到你商店的顾客人数(或者是你的网站的页面访问次数),基于某些确定的特征 X X X ,例如商店的促销、最近的广告、天气、今天周几啊等等。我们已经知道泊松分布(Poisson distribution)通常能适合用来对访客数目进行建模。知道了这个之后,怎么来建立一个模型来解决咱们这个具体问题呢?非常幸运的是,泊松分布是属于指数分布族的一个分部,所以我们可以使用一个广义线性模型(Generalized Linear Model,缩写为 ExpoFamilyGLM)。
对刚刚这类问题如何构建广义线性模型呢?
对于这类问题,我们希望通过一个 X X X 的函数来预测 y \boldsymbol{y} y 的值。为了构建出模型,我们先给出3个假设:
y ∣ x ; θ ∼ E x p o n e n t i a l F a m i l y ( η ) y|x;\theta \sim \rm{ExponentialFamily(\eta)} y ∣ x ; θ ∼ E x p o n e n t i a l F a m i l y ( η ) 。即,给出 x x x 和 η \eta η ,则 y y y 的分布遵循于指数分布。
给出了 x x x 我们的目标是预测 T ( y ) T(y) T ( y ) 的期望值。大多数情况下 T ( y ) = y T(y) = y T ( y ) = y ,也就是说,我们希望通过假设 h h h 输出的 h ( x ) h(x) h ( x ) 能满足 h ( x ) = E [ y ∣ x ] h(x) = E[y|x] h ( x ) = E [ y ∣ x ] 。(统计学知识,有点难)
η \boldsymbol{\eta} η 和 x \boldsymbol{x} x 是线性相关的,
普通最小二乘(Ordinary Least Squares)
普通最小二乘其实是广义线性模型的一个特例,其中 y y y 是连续的,通过 x x x 给出的 y y y 服从高斯分布 N ( 0 , σ 2 ) \mathcal{N}(0,\sigma^2) N ( 0 , σ 2 ) ,经过上面的学习我们有:
h θ = E [ y ∣ x ; θ ] = μ = η = θ T x . \begin{aligned}
h_{\theta}&=E[\boldsymbol{y}|\boldsymbol{x};\boldsymbol{\theta}]\\
&=\boldsymbol{\mu}\\
&=\boldsymbol{\eta}\\
&=\boldsymbol{\theta}^T\boldsymbol{x}.
\end{aligned}
h θ = E [ y ∣ x ; θ ] = μ = η = θ T x .
第一行的等式是基于假设2;第二个等式是基于定理当 y ∣ x ; θ ∼ N ( 0 , σ 2 ) y|x;\theta \sim \mathcal{N}(0,\sigma^2) y ∣ x ; θ ∼ N ( 0 , σ 2 ) ,则 y 的期望就是 μ;第三个等式是基于假设1,以及之前我们此前将高斯分布写成指数族分布的时候推导出来的性质 μ = η \boldsymbol{\mu}=\boldsymbol{\eta}\\ μ = η ;第四个等式就是基于假设3。
Logistic Regression
二分类问题 y ∈ { 0 , 1 } y\in \{0,1\} y ∈ { 0 , 1 } ,可以通过伯努利分布(Bernoulli distribution)来对给定 x x x 的 y y y 进行建模。
伯努利分布,有 ϕ = 1 1 + e − η \phi=\frac{1}{1+e^{-\eta}} ϕ = 1 + e − η 1 ,和在 $ y|x;\theta \sim \rm{Bernoulli}(\phi)$ 下有 E [ y ∣ x ; θ ] = ϕ \quad E[y|x;\theta] = \phi E [ y ∣ x ; θ ] = ϕ 。
则有:
h θ = E [ y ∣ x ; θ ] = ϕ = 1 1 + e − η = 1 1 + e − θ T x . \begin{aligned}
h_{\theta}&=E[\boldsymbol{y}|\boldsymbol{x};\boldsymbol{\theta}]\\
&=\boldsymbol{\phi}\\
&=\frac{1}{1+e^{-\boldsymbol{\eta}}}\\
&=\frac{1}{1+e^{-\boldsymbol{\theta}^T\boldsymbol{x}}}.
\end{aligned}
h θ = E [ y ∣ x ; θ ] = ϕ = 1 + e − η 1 = 1 + e − θ T x 1 .
这就是为什么在 Logistic Function 中我们用 1 1 + e − z \frac{1}{1+e^{-z}} 1 + e − z 1 做假设,即,一旦我们假设给定 x x x 的 y y y 的分布是伯努利分布,那么根据广义线性模型和指数分布族的定义,就会得出这个式子。
注:一个自然参数 η \eta η 的函数 g g g ,g ( η ) = E [ T ( y ) ∣ η ] g(\eta)=E[T(y)|\eta] g ( η ) = E [ T ( y ) ∣ η ] ,这个函数叫做规范响应函数(canonical response function),它的反函数 g − 1 g^{-1} g − 1 叫做规范链接函数(canonical link function)。因此,对于高斯分布来说,它的规范响应函数正好就是识别函数(identify function);而对于伯努利分布来说,它的规范响应函数则是 logistic function。
Softmax Regression
对于多分类问题,有 y ∈ { 1 , 2 , ⋯ , k } y\in \{1,2,\cdots,k\} y ∈ { 1 , 2 , ⋯ , k } ,通过多项式分布(multinomial distribution) 建模。
把多项式推出一个广义线性模型,首先把多项式分布用指数分布族进行描述。
我们给出 k k k 个参数 ϕ 1 , ⋯ , ϕ k \phi_1,\cdots,\phi_k ϕ 1 , ⋯ , ϕ k ,对应各自的输出值的概率,由于 ∑ i = 1 k ϕ i = 1 \sum^k_{i=1}\phi_i=1 ∑ i = 1 k ϕ i = 1 ,所以,有 ϕ k = 1 − ∑ i = 1 k − 1 ϕ i \phi_k =1- \sum^{k-1}_{i=1} \phi_i ϕ k = 1 − ∑ i = 1 k − 1 ϕ i 。注意,ϕ i \phi_i ϕ i 其实是 p ( u = i ; ϕ ) p(u=i;\phi) p ( u = i ; ϕ ) 。现在该出 T ( y ) T(y) T ( y ) :
T ( 1 ) = [ 1 0 0 ⋮ 0 ] , T ( 2 ) = [ 0 1 0 ⋮ 0 ] , T ( 3 ) = [ 0 0 1 ⋮ 0 ] , ⋯ , T ( k − 1 ) = [ 0 0 0 ⋮ 1 ] , T ( k ) = [ 0 0 0 ⋮ 0 ] T(1)=\left[\begin{matrix} 1 \\\\0 \\\\0 \\\\ \vdots\\\\0\end{matrix} \right],
T(2)=\left[\begin{matrix} 0 \\\\1 \\\\0 \\\\ \vdots\\\\0\end{matrix} \right],
T(3)=\left[\begin{matrix} 0 \\\\0 \\\\ 1 \\\\ \vdots\\\\0\end{matrix} \right],\cdots,
T(k-1)=\left[\begin{matrix} 0\\\\0 \\\\0 \\\\ \vdots\\\\1\end{matrix} \right],
T(k)=\left[\begin{matrix} 0\\\\0\\\\0 \\\\ \vdots\\\\0\end{matrix} \right]
T ( 1 ) = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 1 0 0 ⋮ 0 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , T ( 2 ) = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 0 1 0 ⋮ 0 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , T ( 3 ) = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 0 0 1 ⋮ 0 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , ⋯ , T ( k − 1 ) = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 0 0 0 ⋮ 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , T ( k ) = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 0 0 0 ⋮ 0 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤
与之前不同,不再有$ T(y) = y; 然 后 , ;然后, ; 然 后 , T(y)$现在是一个 k – 1 k – 1 k – 1 维的向量,而不是一个实数了。向量 T ( y ) T(y) T ( y ) 中的第 i 个元素写成( T ( y ) ) i (T(y))_i ( T ( y ) ) i 。
给出一个记号:指示函数(indicator function),即 1 { ⋅ } 1\{\cdot\} 1 { ⋅ } 。如果参数为真,则等于1;反之则等于0。
所以我们可以把 T ( y ) T (y) T ( y ) 和 y y y 的关系写成 $ (T(y))_i = 1{y = i}$。
现在把多项式写出指数分布族:
\begin{aligned}
p(y;\theta) &= \phi^{1\\{y=1\\}}\_{1} \phi^{1\\{y=2\\}}\_{2} \cdots \phi^{1\\{y=k\\}}\_{k}\\\\
&= \phi^{1\\{y=1\\}}\_{1} \phi^{1\\{y=2\\}}\_{2} \cdots \phi^{1-\sum^{k-1}\_{i=1}1\\{y=i\\}}\_{k}\\\\
&= \phi^{(T(y))\_1}\_{1} \phi^{(T(y))\_2}\_{2} \cdots \phi^{1-\sum^{k-1}\_{i=1}(T(y))\_i}\_{k}\\\\
&=\rm{exp}\left((T(y))\_1\log{(\phi\_1)}+(T(y))\_2\log{(\phi\_2)} + \cdots + \left(1-\sum^{k-1}\_{i=1}(T(y))\_i\right)\log{(\phi\_k)}\right)\\\\
&=\rm{exp}\left((T(y))\_1\log{(\phi\_1/\phi\_k)}+(T(y))\_2\log{(\phi\_2/\phi\_k)} + \cdots + (T(y))\_{k-1}\log{(\phi\_{k-1}/\phi\_k)} +\log{(\phi\_k)}\right)\\\\
&=b(y)\rm{exp}(\eta^T\it{T}\rm{}\,(y)-a(\eta))
\end{aligned}
其中:
η = [ log ( ϕ 1 / ϕ k ) log ( ϕ 2 / ϕ k ) ⋮ log ( ϕ k − 1 / ϕ k ) ] a ( η ) = − log ϕ i ϕ k b ( y ) = 1 \begin{aligned}
\eta&= \left[\begin{matrix} \log{(\phi_1/\phi_k)}\\\log{(\phi_2/\phi_k)}\\ \vdots\\\log{(\phi_{k-1}/\phi_k)}\end{matrix} \right]\\\\
a(\eta) &= -\log{\frac{\phi_i}{\phi_k}}\\\\
b(y) &= 1
\end{aligned}
η a ( η ) b ( y ) = ⎣ ⎢ ⎢ ⎢ ⎡ log ( ϕ 1 / ϕ k ) log ( ϕ 2 / ϕ k ) ⋮ log ( ϕ k − 1 / ϕ k ) ⎦ ⎥ ⎥ ⎥ ⎤ = − log ϕ k ϕ i = 1
于是对于每一个 η i \eta_i η i 有链接函数:
η i = log ( ϕ i ϕ k ) \eta_i = \log{(\frac{\phi_i}{\phi_k})}
η i = log ( ϕ k ϕ i )
为了简单计算,我们给出定义 η i = log ( ϕ k / ϕ k ) = 0 \eta_i = \log{(\phi_k/\phi_k)} = 0 η i = log ( ϕ k / ϕ k ) = 0 。且对链接函数取反函数然后推导出响应函数,就得到了下面的等式:
e η i = ϕ i ϕ k ϕ k e η i = ϕ i ϕ k ∑ i = 1 k e η i = ∑ i = 1 k ϕ i = 1 \begin{aligned}
e^{\eta_i} & = \frac{\phi_i}{\phi_k}\\\\
\phi_{k}e^{\eta_i} &= \phi_i\\\\
\phi_{k}\sum^k_{i=1}e^{\eta_i}&=\sum^k_{i=1}\phi_i = 1
\end{aligned}
e η i ϕ k e η i ϕ k i = 1 ∑ k e η i = ϕ k ϕ i = ϕ i = i = 1 ∑ k ϕ i = 1
这样得到 ϕ k = 1 / ∑ i = 1 k e η i \phi_k = 1/\sum^k_{i=1}e^{\eta_i} ϕ k = 1 / ∑ i = 1 k e η i ,然后我们我们回代入 e η i = ϕ i ϕ k e^{\eta_i} = \frac{\phi_i}{\phi_k} e η i = ϕ k ϕ i ,
得到相应函数:
ϕ i = e η i ∑ i = 1 k e η i \phi_i = \frac{e^{\eta_i}}{\sum^k_{i=1}e^{\eta_i}}
ϕ i = ∑ i = 1 k e η i e η i
上面这个函数从 η \eta η 映射到了ϕ \phi ϕ ,称为 Softmax函数。通过假设3,我们有了 $\eta_i =\theta_i^Tx $ ,其中的θ 1 , θ 2 , … , θ k − 1 ∈ R n + 1 \theta_1,\theta_2, \dots ,\theta_{k-1} \in \mathbb{R}^{n+1} θ 1 , θ 2 , … , θ k − 1 ∈ R n + 1 就是参数了。我们这里还是定义 θ k = 0 \theta_k=0 θ k = 0 ,这样就有 η k = θ k T x = 0 \eta_k = \theta_k^T x = 0 η k = θ k T x = 0 ,与上文相呼应。
因此,我们有了模型:
p ( y = i ∣ x ; θ ) = ϕ i = e η i ∑ i = 1 k e η i = e θ i T x ∑ i = 1 k e θ i T x \begin{aligned}
p(y=i|x;\theta) &=\phi_i\\
&=\frac{e^{\eta_i}}{\sum^k_{i=1}e^{\eta_i}}\\
&=\frac{e^{\theta_i^Tx}}{\sum^k_{i=1}e^{\theta_i^Tx}}
\end{aligned}
p ( y = i ∣ x ; θ ) = ϕ i = ∑ i = 1 k e η i e η i = ∑ i = 1 k e θ i T x e θ i T x
于是,我们的假设函数是:
h θ ( x ) = E [ T ( y ) ∣ x ; θ ] = E [ 1 { y = 1 } 1 { y = 2 } 1 { y = 3 } x ; θ ⋮ 1 { y = k − 1 } ] = E [ ϕ 1 ϕ 2 ϕ 3 ⋮ ϕ k − 1 ] = E [ e x p θ 1 T x ∑ j = 1 k e x p θ j T x e x p θ 2 T x ∑ j = 1 k e x p θ j T x e x p θ 3 T x ∑ j = 1 k e x p θ j T x ⋮ e x p θ k − 1 T x ∑ j = 1 k e x p θ j T x ] \begin{aligned}
h_\theta(x)&=E[T(y)|x;\theta]\\
&= E \left[\begin{array}{c|c}1\{y = 1\} &\\\\1\{y = 2\}\\1\{y = 3\}& x;\theta\\ \vdots\\1\{y = k-1\}\end{array}\right]\\\\
&=E\left[\begin{array}{c}\phi_1&\\\phi_2\\\phi_3\\ \vdots\\\phi_{k-1}\end{array}\right]\\
&=E\left[\begin{array}{c}
\frac{\rm{exp}^{\theta_1^Tx}}{\sum^k_{j=1}\rm{exp}^{\theta_j^Tx}}&\\
\frac{\rm{exp}^{\theta_2^Tx}}{\sum^k_{j=1}\rm{exp}^{\theta_j^Tx}}\\
\frac{\rm{exp}^{\theta_3^Tx}}{\sum^k_{j=1}\rm{exp}^{\theta_j^Tx}}\\\vdots\\
\frac{\rm{exp}^{\theta_{k-1}^Tx}}{\sum^k_{j=1}\rm{exp}^{\theta_j^Tx}}\end{array}\right]
\end{aligned}
h θ ( x ) = E [ T ( y ) ∣ x ; θ ] = E ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 1 { y = 1 } 1 { y = 2 } 1 { y = 3 } ⋮ 1 { y = k − 1 } x ; θ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ = E ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ϕ 1 ϕ 2 ϕ 3 ⋮ ϕ k − 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ = E ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ∑ j = 1 k e x p θ j T x e x p θ 1 T x ∑ j = 1 k e x p θ j T x e x p θ 2 T x ∑ j = 1 k e x p θ j T x e x p θ 3 T x ⋮ ∑ j = 1 k e x p θ j T x e x p θ k − 1 T x ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤
然后,对于一个训练集来说,我们为了求得 θ \boldsymbol{\theta} θ ,写出似然函数的对数:
ℓ ( θ ) = ∑ i = 1 m log p ( y ( i ) ∣ x ( i ) ; θ ) = ∑ i = 1 m log ∏ l = 1 k ( e x p θ l T x ( i ) ∑ j = 1 k e x p θ j T x ( i ) ) 1 { y ( i ) = l } \begin{aligned}
\ell(\boldsymbol{\theta}) &= \sum^m_{i=1}\log{p(y^{(i)}|x^{(i)};\theta)}\\
&=\sum^m_{i=1}\log{\prod^k_{l=1}\left(\frac{\rm{exp}^{\theta_{l}^Tx^{(i)}}}{\sum^k_{j=1}\rm{exp}^{\theta_j^Tx^{(i)}}}\right)^{1\{y^{(i)}=l\}}}
\end{aligned}
ℓ ( θ ) = i = 1 ∑ m log p ( y ( i ) ∣ x ( i ) ; θ ) = i = 1 ∑ m log l = 1 ∏ k ( ∑ j = 1 k e x p θ j T x ( i ) e x p θ l T x ( i ) ) 1 { y ( i ) = l }
然后我们可以通过梯度上升法或者牛顿法为求:
θ = arg max ℓ ( θ ) \boldsymbol{\theta} = \mathop{\argmax}_{}{\ell(\theta)}
θ = a r g m a x ℓ ( θ )