广义线性模型

Generalized Linear Model

线性回归模型的基本假设为：
$$\begin{gathered} E(Y|X)=\mu (X) \\ \mu (X)=X^T\beta \end{gathered}$$
由此看来$X$只影响到$Y$的均值，不会影响到$Y$的方差，也就是说$Y$的高斯分布为$N(\mu (X),{\sigma }^2)$，因此$Y$的取值为整个实数域。如果将线性模型扩展到广义的模型，需要确定一个问题，什么样的分布可以使用GLM。

Exponential Distribution Family

指数分布族就是这样一个分布族，可以通过一定的方法，使$X$和$Y$的关系转化为线性模型的关系。这个分布族包含了许多分布：比如Gaussion分布、Bernulli分布、Poisson分布、Gamma分布等。指数分布族的表达形式如下：
$$P_{\theta }(x)=P(\theta ,x)=\exp (\theta \cdot x)f(x)g(\theta )$$
其中$\theta \cdot x$为两个参数的内积$\theta$是这个表征分布的参数。通过变形，可以得到下面的形式：
$$P_{\theta }(x)=\exp \left[\sum _{i=1}^k{\eta }_i(\theta )T_i(\theta )-B(\theta )\right]h(x)$$
我们可以认为$B(\theta )$用来归一化，$h(x)$用来处理离散值或者连续值。
以高斯分布为例：
$$\begin{array}{rl}{P}_{\theta }(x)=& \frac{1}{\sigma \sqrt{2\pi }}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)\\ =& \exp \left(-\left(\frac{x^2}{2{\sigma }^2}+\frac{{\mu }^2}{2{\sigma }^2}-\frac{2x\mu }{2{\sigma }^2}\right)-\log (\sigma \sqrt{2\pi })\right)\\ =& \exp \left(\left[-x^2\cdot \frac{1}{2{\sigma }^2}+x\cdot \frac{\mu }{{\sigma }^2}\right]-\left(\frac{{\mu }^2}{2{\sigma }^2}+\log (\sigma \sqrt{2\pi })\right)\right)\end{array}$$
其中前面一部分就是$\eta (\theta )\cdot T(x)$，后面一部分可以认为$B(x)$（可以将$\sqrt{2\pi }$放到指数的系数上），$h(x)$此时可以认为是1。
$$\begin{array}{r}\exp \left(x\cdot \frac{\mu }{{\sigma }^2}-\frac{{\mu }^2}{2{\sigma }^2}\right)\cdot \left(\frac{\exp (-\frac{x^2}{2{\sigma }^2})}{\sigma \sqrt{2\pi }}\right)\end{array}$$
这个形式就一定程度上可以看出高斯分布是由一个中心伴随它的偏离。
接下来是Bernulli分布：
$$\begin{array}{rl}{P}_{\theta }(x)=& p^x(1-p{)}^{(}1-x)\\ =& \exp \left(x\log (\frac{p}{1-p})+\log (1-p)\right)\end{array}$$
接下来介绍的是Canonical Exponential Family。这种分布族在仅有$\theta$一个未知参数的时候，可以得到很多很好的形状，概率密度形式如下：
$$\begin{array}{r}{f}_{\theta }(y)=\exp \left(\frac{y\cdot \theta -b(\theta )}{\varphi }+c(y,\theta )\right)\end{array}$$
这里我们假设$\varphi$是已知的。以高斯分布作为例子：
$$\begin{array}{rl}{f}_{\theta }(y)=& \exp \left(-\left(\frac{x^2}{2{\sigma }^2}+\frac{{\mu }^2}{2{\sigma }^2}-\frac{2x\mu }{2{\sigma }^2}\right)-\log (\sigma \sqrt{2\pi })\right)\\ =& \exp \left(\frac{-\frac{y^2}{2}-\frac{{\mu }^2}{2}+y\mu }{\varphi }-\log (\sigma \sqrt{2\pi }\right)\\ =& \exp \left(\frac{y\theta -\frac{{\mu }^2}{2}}{\varphi }-\left(\log (\sigma \sqrt{2\pi -\frac{y^2}{2\varphi }}\right)\right)\end{array}$$
这里的$b(\theta )=\frac{y^2}{2}$。
因为我们假设了$\varphi$是已知的，因此在转化的过程中，主需要拼凑出$y\theta$这一项。

Likelihood

需要进行对$\theta$的参数估计时，似然函数是我们首先想到的，假设概率密度的对数似然函数为为$l(\theta )=\log (f_{\theta }(y))$，
$$\begin{array}{r}E\left[\frac{\partial l}{\partial \theta }\right]=0\\ E\left[\frac{{\partial }^{2}l}{\partial {\theta }^2}\right]+E{\left[\frac{\partial l}{\partial \theta }\right]}^2=0\end{array}$$
由此可以得到：
$$\begin{array}{rl}mean=& b^{\prime }(\theta )=\mu \\ Var=& \varphi b^{\prime \prime }(\theta )={\sigma }^2\end{array}$$
这样就产生了指数分布族的参数与其矩的对应关系。
简要证明如下：
$$\begin{array}{rl}\frac{\partial }{\partial \theta }\log (f(\theta ))=& \frac{\frac{\partial f(\theta )}{\partial \theta }}{f(\theta )}\\ E=& {\int }_{-\infty }^{\infty }\frac{\frac{\partial f(\theta )}{\partial \theta }}{f(\theta )}f(\theta )d\theta \\ =& \frac{\partial }{\partial \theta }{\int }_{-\infty }^{\infty }f(\theta )d\theta \\ =& 0\\ l(\theta )=& \frac{y\cdot \theta -b(\theta )}{\varphi }+c(y,\theta )\\ \frac{\partial l}{\partial \theta }=& \frac{y-b^{\prime }(\theta )}{\varphi }\\ E\left[\frac{\partial l}{\partial \theta }\right]=& \frac{E(y)-b^{\prime }(\theta )}{\varphi }\\ \frac{{\partial }^2}{\partial {\theta }^2}\log f(\theta )=& \frac{(\frac{{\partial }^{2}f(\theta ))}{\partial {\theta }^2}f(\theta )-\frac{\partial f(\theta )}{\partial \theta }\frac{\partial f(\theta )}{\partial \theta }}{f_{\theta }^2}\\ E=& \int \frac{{\partial }^2}{\partial {\theta }^2}f(\theta )-\frac{(\frac{\partial f(\theta )}{\partial \theta }{)}^2}{f(\theta )}d\theta \\ \frac{{\partial }^{2}l}{\partial {\theta }^2}+\frac{\partial l}{\partial \theta }=& -\frac{b^{\prime \prime }(\theta )}{\varphi }+{\left(\frac{y-b^{\prime }(\theta )}{\varphi }\right)}^2\\ =& 0\\ b^{\prime \prime }(\theta )=& \frac{var(y)}{\varphi }\end{array}$$
我们接下来讨论的分布都属于Canonical Exponential Distribution Family。

Link Function

这些分布由于性质，$Y$的取值有所限定，比如Bernulli分布的取值为 { 0 , 1 } \{0,1\} {0,1}，Poisson分布的取值为非负整数，这样导致了两者均值为$[0,1],[0,+\infty )$。因此需要一个函数将取值转化到整个实数域$\mathbf{R}$上，这个函数就是Link Function $g(\mu (X))=X^T\beta$。比如指数分布：
$$\begin{array}{rl}\mu =& r\exp \delta t\\ \log (\mu )=& \log (r)+\delta t\end{array}$$
其中$r,\delta$是参数，$t,\mu$是对应的$X,Y$。这样的情况下，我们就可以把它抽象为
$$\begin{array}{r}g(\mu )={\beta }_0+X^T{\beta }_1\end{array}$$
也就是说Link Function的目标就是将取值范围拓展到整个实数轴上就行了，对于一个特定的分布，可以说这样的函数有很多种。以Poisson分布为例$Y|XPoisson(\mu (X))$，Link Function需要把非负实数映射到整个实数域，选择对数函数是一个好的方法：$g(\mu (x))=X^T\beta$。至于Bernulli分布，使用的函数为$g(\mu )=\ln (\frac{\mu }{1-\mu })$（logit），或者是${\Phi }^{-1}(\mu )$（probit）。\par
从函数性质来看，link function$g$需要有如下性质：

1. 连续可导
2. 单调严格增长
3. $Im(g)=\mathcal{R}$
4. $g^{-1}$存在且单调递增

那么我们如何选择这样的Link Function。之前谈到指数分布族的参数$\theta ,b^{\prime }(\theta )=\mu$，而$g(\mu )=X^T\beta$。既然$g$可以有很多种选择，$\theta =g(\mu )=X^T\beta$又何妨。这样就直接建立起Canonical EDF和Link Function的关系：
$$\begin{array}{rl}(b^{\prime })=& g^{-1}\\ b^{\prime -1}=& g\end{array}$$
$\theta ,\mu ,\beta$三者建立起了直接关系：
$$\begin{array}{r}\theta (X)\stackrel{b^{\prime }}{⟷}\mu (X)\stackrel{g(\mu )=X^T\beta }{⟷}\beta \\ \theta (x_i)=b^{\prime -1}(g^{-1}(x_i))=h(x_i^T\beta )=x_i^T\beta \end{array}$$
以Bernulli分布为例：
$$\begin{array}{rl}{p}^y(1-p{)}^{1-y}=& \exp \left(y\log \frac{p}{1-p}+\log (1-p)\right)\\ \theta =& \log \frac{p}{1-p}\\ p=& \frac{e^{\theta }}{1+e^{\theta }}\\ f_{\theta }(y)=& \exp (y\theta -\log (1+e^{\theta }))\\ b^{\prime }(\theta )=& \frac{e^{\theta }}{1+e^{\theta }}=p\\ b^{\prime \prime }(\theta )=& \frac{e^{\theta }}{(1+e^{\theta }{)}^2}=p(1-p)\\ g(\mu )=& b^{\prime -1}(\theta )=\log \frac{p}{1-p}\end{array}$$
因此只要得到Canonical EDF中的$b(\theta )$，就可以得到该分布关于GLM的Link Function。

Optimization

在这里要求的参数和线性模型一样，是$\beta$。这里写出对数似然的公式：
$$\begin{array}{rl}{l}_{l}n(\beta )=& \sum _{i=1}^N\frac{\left(y_{i}h(X_i^T\beta )-b(h(X_i^T\beta ))\right)}{\varphi }+c(y_i,\varphi )\\ =& \frac{{\sum }_{i=1}^N\left(y_i{X}_i^T\beta -b(X_i^T\beta )\right)}{\varphi }\\ =& \sum _{i=1}^N\left(\frac{Y_i{X}_i^T\beta }{\varphi }-\frac{b(X_i^T\beta )}{\varphi }\right)\end{array}$$
这个函数中第一项是线形项，二阶导数为0，不改变函数的凹凸性。从Canonical EDF可以知道$b^{\prime \prime }(\theta )=\frac{Var(y)}{\varphi }>0$，因此这个函数严格凸，使用凸优化的方法就可以很好地求解到唯一最佳参数。通过解析的方式得到闭式解也是可以的。

参考：MIT OCW-Statistics for Applications