广义线性模型

 指数分布族

前面学习了线性回归和logistic回归。我们知道对于\(P(y|x;\theta)\)

若y属于实数，满足高斯分布，得到基于最小二乘法的线性回归,若y取{0,1}，满足伯努利分布，得到Logistic回归。

这两个算法，其实都是广义线性模型的特例。

1) 伯努利分布

http://www.cnblogs.com/mikewolf2002/p/7667944.html

2) 高斯分布

http://www.cnblogs.com/mikewolf2002/p/7669977.html

指数分布族的定义：

     若一类概率分布可以写成如下形式，那么它就属于指数分布族：

\[P(y;\eta)=b(y)exp(\eta^{T}T(y)-a(\eta))\]

     \(\eta\)称作分布的自然参数，通常是一个实数。\(T(y)\)称作充分统计量，通常\(T(y)=y\)，实际上是一个概率分布的充分统计量（统计学知识）。对于给定的\(a，b，T\)三个函数，上式定义了一个以\(\eta\)为参数的概率分布集合，即改变\(\eta\)可以得到不同的概率分布。

https://en.wikipedia.org/w/index.php?title=Exponential_family&oldid=648989632

证明伯努利分布是指数分布族：

伯努利分布式对0,1问题进行建模，它的分布律如下：

\[Ber(\phi ),p(y=1; \varphi ) = \varphi; p(y=0; \varphi ) = 1- \varphi\]

伯努利分布是对0，1问题进行建模的，对于\(Bernoulli（\varphi）,y\in \{0,1\}\)有\(p(y=1;\varphi)=\varphi;p(y=0;\varphi)=1−\varphi\)，下面将其推导成指数分布族形式：

\begin{align*} L(\theta) &=P(y;\varphi)=\varphi^y(1-\varphi)^{1-y}=exp(log(\varphi^y(1-\varphi)^{1-y})) \\ 
&=exp(ylog(\varphi)+(1-y)log(1-\varphi)) \\ &= exp(ylog(\frac{\varphi}{1-\varphi})+log(1-\varphi))\end{align*}

    将其与指数族分布形式对比，可以看出：

\[b(y) = 1\]

\[T(y) = y\]

\[\eta=log\frac{\phi}{1-\phi} => \phi=\frac{1}{1+e^{-\eta}}\]

\[a(\eta)=-log(1-\phi)=log(1+e^{\eta})\]

    表明伯努利分布也是指数分布族的一种。从上述式子可以看到，\(\eta\)的形式与logistic函数（sigmoid）一致，这是因为 logistic模型对问题的前置概率估计其实就是伯努利分布。

高斯分布

    下面对高斯分布进行推导，推导公式如下（为了方便计算，我们将方差\(\sigma^2\)设置为1）：

\[N(\mu,1)=\frac{1}{\sqrt{2\pi}}exp^{-\frac{1}{2}(y-\mu)^2}=\frac{1}{\sqrt{2\pi}}exp^{-\frac{1}{2}y^2+y\mu-\frac{1}{2}\mu^2}=\frac{1}{\sqrt{2\pi}}exp^{-\frac{1}{2}y^2}exp^{y\mu-\frac{1}{2}\mu^2}\]

将上式与指数族分布形式比对，可知：

\[b(y)=\frac{1}{\sqrt{2\pi}}exp^{-\frac{1}{2}y^2}\]

\[T(y)=y\]

\[\eta=\mu\]

\[a(\eta)=\frac{\mu^2}{2}\]

两个典型的指数分布族，伯努利和高斯分布。其实大多数概率分布都可以表示成指数分布族形式，如下所示：

• 伯努利分布（Bernoulli）：对 0、1 问题进行建模；
• 多项式分布（Multinomial）：对 K 个离散结果的事件建模；
• 泊松分布（Poisson）：对计数过程进行建模，比如网站访问量的计数问题，放射性衰变的数目，商店顾客数量等问题；
• 伽马分布（gamma）与指数分布（exponential）：对有间隔的正数进行建模，比如公交车的到站时间问题；
• β 分布：对小数建模；
• Dirichlet 分布：对概率分布进建模；
• Wishart 分布：协方差矩阵的分布；

广义线性模型GLM

选定了一个指数分布族后，怎样来推导出一个GLM呢？

假设：

（1）\(y|x; \theta~ExponentialFamily（\eta）\)，即假设试图预测的变量\(y\)在给定\(x\)，以\(\theta\)作为参数的条件概率，属于以\(\eta\)作为自然参数的指数分布族，例：若要统计网站点击量\(y\)，用泊松分布建模。

（2） 给定\(x\)，目标是求出以\(x\)为条件的\(T(y)\)的期望\(E[T(y)|x]\)，即让学习算法输出\(h(x) = E[T(y)|x]\)

（3）\(\eta=\theta^Tx\),即自然参数和输入特征\(x\)之间线性相关，关系由θ决定。仅当\(\eta\)是实数时才有意义,若\(\eta\)是一个向量，则\(\eta_i=\theta_i^Tx\)

推导伯努利分布的GLM：

\(y|x; \theta~ExponentialFamily（\eta）\)，伯努利分布属于指数分布族,对给定的\(x,\theta\)，学习算法进行一次预测的输出：

\[h_\theta(x)=E[y|x;\theta]=\phi=\frac{1}{1+e^{-\eta}}=\frac{1}{(1+e^{-\theta^Tx})}\]

得到logistic回归算法。

正则响应函数：\(g(\eta)=E[y;\eta]\),将自然参数\(\eta\)和\(y\)的期望联系起来

正则关联函数:\(g^{-1}\)

推导多项式分布的GLM：

多项式分布是在\(k\)个可能取值上的分布，即\(y\in {1,..,k}\),如将收到的邮件分成\(k\)类，诊断某病人为\(k\)种病中的一种等问题。

（1）将多项式分布写成指数分布族的形式：

设多项式分布的参数：\(\phi_1,\phi_2,…,\phi_k\) ，且\(\sum\limits_{i=1}^k\phi_i=1\) ，\(\phi_i\)表示第\(i\)个取值的概率分布，最后一个参数可以由前\(k-1\)个推导出，所以只将前\(k-1\)个\(\phi_1,\phi_2,…,\phi_{k-1}\)视为参数。

多项式分布是少数几个\(T(y)!=y\)的分布，\(T(1) \sim T(k)\)都定义成一个\(k-1\)维的列向量，表示为：

\[T(1)=\begin{bmatrix}
1\\
0\\
0\\
\vdots\\
0\end{bmatrix}T(2)=\begin{bmatrix}
0\\
1\\
0\\
\vdots\\
0\end{bmatrix}T(3)=\begin{bmatrix}
0\\
0\\
1\\
\vdots\\
0\end{bmatrix}T(k-1)=\begin{bmatrix}
0\\
0\\
0\\
\vdots\\
1\end{bmatrix}T(k)=\begin{bmatrix}
0\\
0\\
0\\
\vdots\\
0\end{bmatrix}\]

这样定义\(T(y)\)是为了将多项式分布写成指数分布族形式。

*定义符号：指示函数，\(1\{.\}\)

1{True} = 1, 1{False} = 0，即大括号内命题为真，值为1,否则为0。

例：1{2=3} = 0, 1{1+1=2} = 1

用\(T(y)_i\)表示\(T(y)\)的第\(i\)个元素，则\(T(y)_i=1\{y=i\}\)

根据参数\(\phi\)的意义（\(\phi_i\)表示第\(i\)个取值的概率分布），可推出：

\begin{align*} p(y;\phi) &=\phi_1^{1\{y=1\}}\phi_2^{1\{y=2\}}...\phi_k^{1\{y=k\}} \\ 
&=\phi_1^{1\{y=1\}}\phi_2^{1\{y=2\}}...\phi_k^{1-\sum\limits_{i=1}^{k-1}1\{y=i\}} \\ &=\phi_1^{(T(y))_1}\phi_2^{(T(y))_2}...\phi_k^{1-\sum\limits_{i=1}^{k-1}(T(y))_i} \\ &=>exp(T(y))_1log(\phi_1)+ exp(T(y))_2log(\phi_2)+...+(1-\sum\limits_{i=1}^{k-1}(T(y))_i)log(\phi_k) \\ &=exp(T(y))_1log(\phi_1/\phi_k)+ exp(T(y))_2log(\phi_2/\phi_k)+...+exp(T(y))_{k-1}log(\phi_{k-1}/\phi_k)+log(\phi_k)) \\ &=b(y)exp(\eta^{T}T(y)-a(\eta))\end{align*}

可得：

\[\eta=\begin{bmatrix}
log(\phi_1/\phi_k)\\
log(\phi_2/\phi_k)\\
\vdots\\
log(\phi_{k-1}/\phi_k)
\end{bmatrix}\]

\[a(\eta)=-log(\phi_k)\]

\[b(y)=1\]

证明多项式分布是指数分布族。

再用\(\eta\)表示\(\phi\)：

\[\eta_i=log\frac{\phi_i}{\phi_k} \\
e^{\eta_i}=\frac{\phi_i}{\phi_k} \\
\phi_ke^{\eta_i}=\phi_i \\
\phi_k\sum\limits_{i=1}^ke^{\eta_i}=\sum\limits_{i=1}^k\phi_i=1 \\
\phi_i=\frac{e^{\eta_i}}{\sum\limits_{j=1}^ke^{\eta_j}}\]

（2）根据上述假设（3）中自然参数和输入\(x\)的线性关系，可求得：

\begin{align*} p(y=i|x;\theta) &=\phi_i \\ 
  &=\frac{e^{\eta_i}}{\sum\limits_{j=1}^ke^{\eta_j}} \\ &=\frac{e^{\theta_i^Tx}}{\sum\limits_{j=1}^ke^{\theta_j^Tx}} \\ \end{align*}

（3）根据上述假设（2）中的输出\(h(x) = E[T(y)|x]\)，可求得：

\begin{align*} h_\theta(x) &=E[T(y)|x;\theta] \\ 
  &=\begin{bmatrix}
1\{y=1\}|x;\theta\\
1\{y=2\}|x;\theta\\
\vdots\\
1\{y=k-1\}|x;\theta
\end{bmatrix} \\ &= \begin{bmatrix}
\phi_1\\
\phi_2\\
\vdots\\
\phi_{k-1}
\end{bmatrix}\\ &=\begin{bmatrix}
\frac{e^{\theta_1^Tx}}{\sum\limits_{j=1}^ke^{\theta_j^Tx}} \\
\frac{e^{\theta_2^Tx}}{\sum\limits_{j=1}^ke^{\theta_j^Tx}} \\
\vdots\\
\frac{e^{\theta_{k-1}^Tx}}{\sum\limits_{j=1}^ke^{\theta_j^Tx}}
\end{bmatrix} \end{align*}

称这种回归算法为softmax回归，是logistic回归的推广。

Softmax回归的训练方法和logistic相似，通过极大似然估计找到参数\(\theta\),其对数似然性为：

\begin{align*} l(\theta) &=\sum\limits_{i=1}^{m}logp(y^{(i)}|x^{(i)};\theta); \\ 
  &=\sum\limits_{i=1}^{m}log\prod\limits_{l=1}^k\bigg(\frac{e^{\theta_l^Tx^{(i)}}}{\sum\limits_{j=1}^ke^{\theta_j^Tx^{(i)}}}\bigg)^{1\{y^{(i)}=l\}} = \sum\limits_{i=1}^m\bigg(\sum\limits_{l=1}^k\bigg( 1\{y^{(i)}=l\}\log\bigg( \frac{e^{\theta_l^Tx^{(i)}}}{\sum\limits_{j=1}^ke^{\theta_j^Tx^{(i)}} }\bigg)\bigg)\bigg)\end{align*}

其中\(\theta\)是一个\(k\times(n+1)\)的矩阵，代表这\(k\)个类的所有训练参数，每个类的参数是一个\(n+1\)维的向量。

跟Logistic回归一样，softmax也可以用梯度下降法或者牛顿迭代法求解，但在softmax回归中直接用上述对数似然函数求导是不能更新参数的，因为它存在冗余的参数，通常用牛顿方法中的Hessian矩阵也不可逆，是一个非凸函数，那么可以通过添加一个权重衰减项来修改代价函数，使得代价函数是凸函数，并且得到的Hessian矩阵可逆

再通过梯度上升或牛顿方法找对数似然性的最大值，即可求出参数\(\theta\)。

详细的公式推导见：http://ufldl.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92