GLM 广义线性模型

George Box said: “All models are wrong, some are useful”

1. 始于 Linear Model

作为 GLM 的基础，本节 review 经典的 Linear Regression，并阐述一些基础 term。
我们线性回归的基本如下述公式，本质上是想通过观察 xx，然后以一个简单的线性函数 h(x)h(x) 来预测 yy：

y=h(x)=wTx

y=h(x)=w^Tx

1.1 dependent variable yy

这是我们的预测目标，也称 response variable。这里有一个容易混淆的点，实际上 yy 可以表达三种含义（建模用的是分布，观察到的是采样，预测的是期望）：

distribution；抽象地讨论 response variable 时，我们实际上关注对于给定数据和参数时， y| x,wy|\ x,w 服从的分布。Linear Regression 的 yy 服从高斯分布，具体取值是实数，但这里我们关注的是分布。
observed outcome；我们的 label，有时用 tt 区分表示；这是真正观察到的结果，只是一个值。
expected outcome；y=E[y|x]=h(x)y=\mathbf{E}[y|x]=h(x) 表示模型的预测；注意 yy 实际上服从一个分布，但预测结果是整个分布的均值 μ\mu，只是一个值。

1.2 independent variable xx

这是我们的特征，可以包含很多维度，一个特征也称为一个 predictor。

1.3 hypothesis h(x)h(x)

线性模型的假设非常简单，即 h(x)=wTxh(x) = w^Tx inner product of weight vector and feature vector，被称为 linear predictor。这是就是线性模型，GLM 也是基与此的推广。
深入来看，各个维度特征（predictor） xjx_j 通过系数 wjw_j 线性加和，这一过程将信息进行了整合；而不同的 weight（coefficient）反映了相关特征不同的贡献程度。

2. 推广到 Generalized Linear Model

2.1 Motive & Definition

线性模型有着非常强的局限，即 response variable yy 必须服从高斯分布；主要的局限是拟合目标 yy 的 scale 是一个实数 (−∞,+∞)(-\infty,+\infty)。具体来说有俩个问题：

yy 的取值范围和一些常见问题不匹配。例如 count（游客人数统计恒为正）以及 binary（某个二分类问题）
yy 的方差是常数 constant。有些问题上方差可能依赖 yy 的均值，例如我预测目标值越大方也越大（预测越不精确）

所以这时我们使用 Generalized Linear Model 来克服这俩个问题。
一句话定义 GLM 即（from wiki）：
In statistics, the generalized linear model (GLM) is a flexible generalization of ordinary linear regression that allows for response variables that have error distribution models other than a normal distribution.
详细来说，我们可以把 GLM 分解为 Random Component、System Component 和 Link Function 三个部分。

2.2 Random Component

An exponential family model for the response

这里是指 response variable 必须服从某一 exponential family distribution 指数族分布，即 y|x,w∼ExponentialFamily(η)y|x,w \thicksim ExponentialFamily(\eta)，η\eta 指 exponential family 的 natural parameter 自然参数。
例如 linear regression 服从 Gaussian 高斯分布，logistic regression 服从 Bernoulli 伯努利分布。指数族还有很多分布如多项分布、拉普拉斯分布、泊松分布等等。
另外，这也可以被称为 Error Structure : error distribution model for the response。对于 Gaussian 的 residual 残差 ϵ=y−h(x)\epsilon=y-h(x) 服从高斯分布 N(0,σ)N(0,\sigma) 是很直观可；但是，例如 Bernoulli 则没有直接的 error term。可以构造其他的 residual 服从 Binomial，但是较抽象，本文不对此角度展开。

2.3 Systematic Component

linear predictor

广义线性模型 GLM 本质上还是线性模型，我们推广的只是 response variable yy 的分布，模型最终学习的目标还是 linear predictor wTxw^Tx 中的 weight vector。
注意，GLM 的一个较强的假设是 η=wTx\eta=w^Tx，即 yy 相关的 exponential family 的 natural parameter η\eta 等于 linear predictor。这个假设倾向于是一种 design choice（from Andrew），不过这种 choice 是有合理性的，至少 η\eta 和 linear predictor 的 scale 是一致的。

2.4 Link Function

A link function connects the mean of the response to the linear predictor

通过上述的 Random Component 和 Systematic Component，我们已经把 yy 和 wTxw^Tx 统一到了 exponential family distribution 中，最终的一步就是通过 link function 建立俩者联系。对任意 exponential family distribution，都存在 link function g(μ)=ηg(\mu)=\eta，μ\mu 是分布的均值而 η\eta 是 natural parameter；例如 Gaussian 的 link function 是 identity（g(μ)=ηg(\mu)=\eta），Bernoulli 的 link function 是 logit（g(μ)=lnμ1−μ=ηg(\mu)=ln\displaystyle\frac{\mu}{1-\mu}=\eta）。
link function 建立了response variable 分布均值（实际就是我们的预测目标）和 linear predictor 的关系（准确来说，这只在 T(y)=yT(y)=y 条件下成立，但大多数情况都条件都成立，这里不展开说明）。实际上 link function 把原始 yy 的 scale 转换统一到了 linear predictor 的 scale 上。另外，不同分布的 link function 可以通过原始分布公式变换到指数组分布形式来直接推出，之后本文第4章会有详细讲解。
最后要强调的是，link function 的反函数 g−1(η)=μg^{-1}(\eta)=\mu 称为响应函数 response function。响应函数把 linear predictor 直接映射到了预测目标 yy，较常用的响应函数例如 logistic/sigmoid、softmax（都是 logit 的反函数）。

2.5 Contrast between LM & GLM

linear predictor η=wTx\eta=w^Tx

Linear Regression

response variable y∼N(η,σ2e)y \thicksim N(\eta,\sigma_e^2)
link function η=g(μ)=μ\eta=g(\mu)=\mu, called identity
prediction h(x)=E[y|x,w]=μ=g−1(η)=μh(x) = \mathbf{E}[y|x,w]=\mu=g^{-1}(\eta)=\mu

Generalized Linear Model

response variable y∼exponential familyy \thicksim exponential\ family
link function g(μ)g(\mu), eg. logit for Bernoulli
prediction h(x)=E[y|x,w]=μ=g−1(η)h(x) = \mathbf{E}[y|x,w]=\mu=g^{-1}(\eta), eg.logistic for Bernoulli

这里再次强调了他们 linear predictor 的部分是一致的；不过对 response variable 服从分布的假设不一致。Gaussian 的 response function 是 g−1(η)=μg^{-1}(\eta)=\mu；而 exponential family 根据具体假设的分布，使用相应的 response function （例如 Bernoulli 是 sigmoid）。

额外强调一个点，无论是 LM 还是 GLM，我们对不同数据 xx 得到的其实是不同的 response variable 的分布，所以不同分布的 μ\mu 值不同，进而我们预测的结果不同。虽然每一条数据只是预测了一个值，但其实对应的是一个分布。并且数据相似的话对应的分布也相似，那么预测结果也相似。

3. 实例：Linear Regression 到 Logistic Regression

GLM 中我们最常用到的是 Logistic Regression；即假设 yy 服从 Bernoulli 伯努利分布，这里详细展开阐述。

3.1 以 GLM 来看 Logistic Regression

以下直接列出相关的 GLM 概念对 LR 的解释：

An exponential family (random component) yi∼Bern(μi)y_i\thicksim Bern(\mu_i)
- linearl predictor (system component)：ηi=∑Jjwjxij\eta_i=\sum_j^J w_jx_{ij}，这个部分 GLM 都有且都一致
- link function：η=g(μ)=lnμ1−μ\eta=g(\mu)=ln\displaystyle\frac{\mu}{1-\mu}，这个函数是 log-odds，又称 logit 函数
- response function：μ=g−1(η)=11+e−η\mu=g^{-1}(\eta)=\displaystyle\frac{1}{1+e^{-\eta}}，称为 logistic 或 sigmoid
- prediction：h(xi)=E[yi|w,xi]=sigmoid(ηi)h(x_i)=\mathbf{E}[y_i|w,x_i]=sigmoid(\eta_i)
- loss function : E=−ylnh(x)−(1−y)ln(1−h(x))E=-ylnh(x)-(1-y)ln(1-h(x)), 和 Linear Model 一样由 MLE 推导得到，3.3 中会给出详细推导过程
- Scale Insight 可以给我们更多 link function 的理解：
  - binary outcome：0 or 1，对应 observed outcome 即 label
  - probability：[0,1][0,1]，对应 expected outcome，即 yy 的分布的 μ\mu
  - odds：(0,∞)(0,\infty)，概率和几率可相互转换（oa=pa1−pa, pa=oa1+oao_a=\displaystyle\frac{p_a}{1-p_a},\ p_a=\displaystyle\frac{o_a}{1+o_a}）；即是发生概率与不发生概率的比，赌博中常见这一概念，和赔率相关
  - log-odds/logit：(−∞,+∞)(-\infty,+\infty)，即 logpa1−palog\displaystyle\frac{p_a}{1-p_a}
  所以 log-odds 在 scale 上和 linear predictor 匹配。对于 Logistic Regression，我们通过 link function logit 建立了与简单线性模型的关联。link function 就是把任意 exponential family distribution 的均值 μ\mu 映射到线性尺度上：
  
  η=g(μ)=lnμ1−μη=lnE[yi|w,xi]1−E[yi|w,xi]=wTxi=∑jJwjxij
  
  \eta=g(\mu)=ln\displaystyle\frac{\mu}{1-\mu}\\\eta=ln\frac{\mathbf{E}[y_i|w,x_i]}{1-\mathbf{E}[y_i|w,x_i]}=w^Tx_i=\sum_j^J w_jx_{ij}
  即是，在 Logistic Regression 模型中，输出 y=1y=1 的 log-odds 对数几率是关于输入 xx 的线性函数。
  
  3.2 Logistic/Sigmoid 函数
  
  此函数是 Bernoulli 在 GLM 的 response function，即 link function 的反函数。（强调反函数的时候可以使用 g−1(x)g^{-1}(x) 标识；当然直接使用 g(x)g(x) 也没问题，Andrew 课上使用的是后者）
  此函数在神经网络中也有用作 activate function 来引入非线性（当然现在更常用 rectified linear unit, RELU），其拥有一些特点：
  - bounded 有界的，取值 (0,1)(0,1)，有概率意义
  - monotonic 单调递增函数
  - differentiable 可微，且求导简单
  - 函数呈 S 曲线，中央区域敏感，两侧抑制（有俩个优点：首先是符合神经激活原理；其次是现实问题会有类似的非线性关系，当自变量很大或很小时对因变量影响很小，但是在某个阈值范围上时影响很大）
  这里展开证明一下求导简单，即 g′(z)=g(z)(1−g(z))g'(z)=g(z)(1-g(z))：
  
  sigmoid=g(z)g′(z)=11+e−z=−(e−z∗−1)(1+e−z)2=e−z(1+e−z)2=(1+e−z)−1(1+e−z)2=11+e−z−1(1+e−z)2=g(z)−g(z)2=g(z)(1−g(z))
  
  \begin{aligned} sigmoid=g(z)&=\frac{1}{1+e^{-z}}\\ g'(z)&=\frac{-(e^{-z}*-1)}{(1+e^{-z})^2}=\frac{e^{-z}}{(1+e^{-z})^2}\\ &=\frac{(1+e^{-z})-1}{(1+e^{-z})^2}=\frac{1}{1+e^{-z}}-\frac{1}{(1+e^{-z})^2}\\ &=g(z)-g(z)^2=g(z)(1-g(z)) \end{aligned}
  
  另外，这里先证明一个有趣的特性 g(−z)=1−g(z)g(-z)=1-g(z)（为后续章节证明流程铺垫~）
  
  1−g(z)=1−11+e−z=e−z1+e−z=11+ez=g(−z)
  
  1-g(z)=1-\frac{1}{1+e^{-z}}=\frac{e^{-z}}{1+e^{-z}}=\frac{1}{1+e^{z}}=g(-z)
  
  logistic/sigmoid 函数还有一些其他的解释，例如生态学模型。
  
  3.3 LR 的 loss function
  
  3.3.1 标准套路
  
  本节详细完成 logistic regression 剩下部分的阐述。假设我们观察到了数据集 (xi,ti)(x_i, t_i) 包含 NN 个样本，xix_i 表示某一样本的特征向量，ti=0 or 1t_i=0\ or\ 1 表示这一样本的真实类别（observed outcome）。另外定义 yi=h(xi)y_i=h(x_i) 为模型的输出结果（expected outcome）。我们假设 tit_i 服从 Bernoulli 伯努利分布（Bern(x|μ)=μx(1−μ)1−xBern(x|\mu)=\mu^{x}(1-\mu)^{1-x}），则可以写出似然函数如下：
  
  L(D|w)=p(t|w)=∏i=1Nytii(1−yi)1−ti=∏i=1Nh(xi)ti(1−h(xi))1−ti
  
  L(D|w)=p(t|w)=\prod_{i=1}^Ny_i^{t_i}(1-y_i)^{1-t_i}=\prod_{i=1}^Nh(x_i)^{t_i}(1-h(x_i))^{1-t_i}
  对似然函数取对数化简如下：
  
  l=lnp(t|w)=∑i=1N{tilnh(xi)+(1−ti)ln(1−h(xi))}
  
  l=lnp(t|w)=\sum_{i=1}^N\{t_ilnh(x_i)+(1-t_i)ln(1-h(x_i))\}
  所以 loss fucntion 可以写为（最大化似然转最小化cost）：
  
  J=−l=∑i=1N{−tilnh(xi)−(1−ti)ln(1−h(xi))}
  
  J=-l=\sum_{i=1}^N\{-t_ilnh(x_i)-(1-t_i)ln(1-h(x_i))\} 这是基于 binomial 二项分布推导出来的 loss，所以又被称为 negative binomial log-likelihood。
  
  其中 h(x)h(x) 依照 GLM 理论，依据 response function 展开为：
  
  h(x)=g−1(wTx)=11+e−wTx
  
  h(x)=g^{-1}(w^Tx)=\frac{1}{1+e^{-w^Tx}}
  为 loss function JJ 带入 h(xi)h(x_i) 对 weight vector ww 求导（注意 Chain Rule）：
  
  ∂ Ji∂ w=∂ Ji∂ h(xi)∂ h(xi)∂ w=(−ti1h(xi)−(1−ti)11−h(xi)∗(−1))∗(h(xi)(1−h(xi))xi)=−ti1h(xi)h(xi)(1−h(xi))xi+(1−ti)11−h(xi)h(xi)(1−h(xi))xi=−ti(1−h(xi))xi+(1−ti)h(xi)xi=−tixi+tih(xi)xi+h(xi)xi−tih(xi)xi=(h(xi)−ti)xi
  
  \begin{aligned} \frac{\partial\ J_i}{\partial\ w}&=\frac{\partial\ J_i}{\partial\ h(x_i)}\frac{\partial\ h(x_i)}{\partial\ w}\\ &=(-t_i\frac{1}{h(x_i)}-(1-t_i)\frac{1}{1-h(x_i)}*(-1))*(h(x_i)(1-h(x_i))x_i)\\ &=-t_i\frac{1}{h(x_i)}h(x_i)(1-h(x_i))x_i+(1-t_i)\frac{1}{1-h(x_i)}h(x_i)(1-h(x_i))x_i\\ &=-t_i(1-h(x_i))x_i+(1-t_i)h(x_i)x_i\\ &=-t_ix_i+t_ih(x_i)x_i+h(x_i)x_i-t_ih(x_i)x_i\\ &=(h(x_i)-t_i)x_i \end{aligned}
  所以 Logistic Regression 与 Linear Regression 更新公式（Batch Gradient Descent）是类似的，唯一不同点在于不同的 response function h(x)h(x) ：
  
  wnew=w−λ∂ J∂ w=w−λ∑i=1N(h(xi)−ti)xi=w−λ∑i=1N(yi−ti)xi
  
  w^{new}=w-\lambda\frac{\partial\ J}{\partial\ w}=w-\lambda\sum_{i=1}^{N}(h(x_i)-t_i)x_i=w-\lambda\sum_{i=1}^{N}(y_i-t_i)x_i
  
  3.3.2 其他套路
  
  这里专门强调一下，有一些论文里面会使用其他形式的等价的 negative binomial log-likelihood：
  
  J=∑i=1Nln(1+exp{−tizi})
  
  J=\sum_{i=1}^Nln(1+\exp\{-t_iz_i\})
  这个形式要简洁很多，但实际和上小节中 loss function J=∑Ni=1{−tilnh(xi)−(1−ti)ln(1−h(xi))}J=\sum_{i=1}^N\{-t_ilnh(x_i)-(1-t_i)ln(1-h(x_i))\} 完全等价的，唯一的区别是：这里使用的 label ti=−1/+1t_i=-1/+1 而非上一小节的 ti=0/1t_i=0/1。这种表达有一个好处， ziz_i 是 linear predictor，上小节有 zi=wTxiz_i=w^Tx_i，但是这里更 general，可以用其他方式定义来定义这个 linear predictor，例如 GBDT 中的函数加和，其实这个 loss 就是从 GBDT 论文上看到的~以下将证明这一套路与标准讨论等价：
  
  首先，把原 loss 展开 h(xi)=g−1(wTxi)=g(zi)h(x_i)=g^{-1}(w^Tx_i)=g(z_i)（把 inverse 符号去掉了，这里不需要再强调了），且根据 sigmoid 的 g(−z)=1−g(z)g(-z)=1-g(z)，有：
  
  J=∑i=1N{−tilnh(xi)−(1−ti)ln(1−h(xi))}=∑i=1N{−tiln(g(zi))−(1−ti)ln(g(−zi))}
  
  J=\sum_{i=1}^N\{-t_ilnh(x_i)-(1-t_i)ln(1-h(x_i))\}\\ =\sum_{i=1}^N\{-t_iln(g(z_i))-(1-t_i)ln(g(-z_i))\}
  对 label tit_i 分情况展开：
  
  Ji={−ln(g(zi))−ln(g(−zi))if ti=1if ti=0
  
  J_i=\begin{cases} -ln(g(z_i))& if\ t_i=1\\ -ln(g(-z_i))& if\ t_i=0 \end{cases}
  切换到 ti=−1/+1t_i=-1/+1（这是唯一的不同哦）：
  
  Ji={−ln(g(tizi))−ln(g(tizi))if ti=+1if ti=−1
  
  J_i=\begin{cases} -ln(g(t_iz_i))& if\ t_i=+1\\ -ln(g(t_iz_i))& if\ t_i=-1 \end{cases}
  显然可以对原 loss 整合并展开为：
  
  J=∑i=1N{−ln(g(tizi))}=∑i=1N{−ln(11+exp{−tizi})}=∑i=1N{ln(1+exp{−tizi})}
  
  J=\sum_{i=1}^N\{-ln(g(t_iz_i))\}=\sum_{i=1}^N\{-ln(\frac{1}{1+\exp\{-t_iz_i\}})\}=\sum_{i=1}^N\{ln(1+\exp\{-t_iz_i\})\} 证毕。
  
  loss 的 gradient 当然也是完全等价的，这里简单展开部分结果：
  
  ∂Ji∂z=11+exp{−tizi}exp{−tizi}(−ti)=−ti1+exp{tizi}
  
  \frac{\partial J_i}{\partial z}=\frac{1}{1+\exp\{-t_iz_i\}}\exp\{-t_iz_i\}(-t_i)=\frac{-t_i}{1+\exp\{t_iz_i\}}
  
  4. Exponential Family
  
  这里补充一下 exponential family 的定义，并且给出 Bernoulli 以及 Categorical 对应的 link function 和 response function。
  
  4.1 Definition
  
  The exponential family of distribution over xx given η\eta is （h,g,uh, g, u 都是确定的函数，基于 η\eta 或 xx 确定的，所以这是只以 η\eta 为参数的分布）:
  
  p(x|η)=h(x)g(η)exp{ηTu(x)}
  
  p(x|\eta)=h(x)g(\eta)exp\{\eta^Tu(x)\}
  - η\eta natural parameter 自然参数，这是决定分布的具体参数
  - u(x)u(x) sufficient statistics 充分统计量，通常有 u(x)=xu(x)=x
  - g(η)g(\eta) 是分布正规化系数 coefficient，即确保概率和为1，满足 g(η)∫h(x)exp{ηTu(x)}dx=1g(\eta)\int h(x)exp\{\eta^Tu(x)\}dx = 1
  常见的 exponential family 有 Bernoulli Distribution, Binomial Poisson Distribution, Negative Binomial Distribution, Categorical Distribution, Multinomial Distribution, Beta Distribution, Dirichlet Distribution, Laplace Distribution, Gamma Distribution, Normal Distribution 等等，所以说 GLM 极大地拓展了 LM 的使用范围。
  
  4.2 Bernoulli Distribution
  
  伯努利分布的分布律（相对于 continuous variable 的概率密度函数）如下：
  
  p(x|μ)=Bern(x|μ)=μx(1−μ)1−x
  
  p(x|\mu)=Bern(x|\mu)=\mu^x(1-\mu)^{1-x}
  我们将其变换到标准 exponential family 形式：
  
  p(x|μ)=exp{ln{μx(1−μ)1−x}}=exp{xlnμ+(1−x)ln(1−μ)}=exp{xlnμ−xln(1−μ)+ln(1−μ)}=(1−μ)exp{xlnμ1−μ}
  
  \begin{aligned} p(x|\mu) &=exp\{ln\{\mu^x(1-\mu)^{1-x}\}\}\\ &=exp\{xln\mu+(1-x)ln(1-\mu)\}\\ &=exp\{xln\mu-xln(1-\mu)+ln(1-\mu)\}\\ &=(1-\mu)exp\{xln\frac{\mu}{1-\mu}\} \end{aligned}
  得到 natural parameter，也正是 GLM 使用的 link function，logit：
  
  η=lnμ1−μ
  
  \eta=ln\frac{\mu}{1-\mu}
  这里推导一下其 inverse function，即 GLM 使用的 response function：
  
  e−η=exp{−lnμ1−μ}=1−μμ1+e−η=1+1−μμ=1μ11+e−η=μ
  
  e^{-\eta}=exp\{-ln\frac{\mu}{1-\mu}\}=\frac{1-\mu}{\mu}\\ 1+e^{-\eta}=1+\frac{1-\mu}{\mu}=\frac{1}{\mu}\\ \frac{1}{1+e^{-\eta}}=\mu
  这也就是 logistic/sigmoid 函数。剩余其他部分：
  
  u(x)=x, h(x)=1,g(η)=sigmoid(−η)=11+eη=11+μ1−μ=(1−μ)
  
  u(x)=x,\ h(x)=1,\\g(\eta)=sigmoid(-\eta)=\frac{1}{1+e^{\eta}}=\frac{1}{1+\frac{\mu}{1-\mu}}=(1-\mu)
  
  4.3 Categorical Distribution
  
  作为 Bernoulli 的推广，Categorical 的 link function 和 response function 与前者非常相似。其 response function 是 softmax，所以 logistic/sigmoid 和 softmax 也是一个推广关系。
  
  这里注意一点，Categorical Distribution 即是单次的 Multiple Distribution，后者更常见。（而 Bernoulli Distribution 是单次的 Binomial Distribution）
  
  以下介绍推导过程，分类分布的分布律，以及 exponential family 形式如下：
  
  p(x|μ)=∏k=1Mμxkk=exp{∑k=1Mxklnμk}
  
  p(x|\mu)=\prod_{k=1}^M\mu_k^{x_k}=exp\{\sum_{k=1}^Mx_kln\mu_k\}
  上述表达缺少了一个约束： ∑k=1Mμk=1\displaystyle\sum_{k=1}^M\mu_k=1，通常会改写分布形式来消除这个约束，即我们只确认 M−1M-1 个 μk\mu_k，剩下 1 个 μk\mu_k 是一个确定的值。当然，我们其实还会有隐含的约束 0≤μk≤10\leq\mu_k\leq1 和 ∑k=1M−1μk≤1\displaystyle\sum_{k=1}^{M-1}\mu_k\leq1，这个 Bernoulli 也有。
  下面是包含约束的改写过程：
  
  exp{∑k=1Mxklnμk}=exp{∑k=1M−1xklnμk+(1−∑k=1M−1xk)ln(1−∑k=1M−1μk)}=exp{∑k=1M−1xklnμk−∑k=1M−1xkln(1−∑k=1M−1μk)+ln(1−∑k=1M−1μk)}=exp{∑k=1M−1xklnμk1−∑M−1k=1μk+ln(1−∑k=1M−1μk)}=(1−∑k=1M−1μk)exp{∑k=1M−1xklnμk1−∑M−1k=1μk}
  
  \begin{aligned} exp\{\sum_{k=1}^Mx_kln\mu_k\} &=exp\{\sum_{k=1}^{M-1}x_kln\mu_k+(1-\sum_{k=1}^{M-1}x_k)ln(1-\sum_{k=1}^{M-1}\mu_k)\}\\ &=exp\{\sum_{k=1}^{M-1}x_kln\mu_k-\sum_{k=1}^{M-1}x_kln(1-\sum_{k=1}^{M-1}\mu_k)+ln(1-\sum_{k=1}^{M-1}\mu_k)\}\\ &=exp\{\sum_{k=1}^{M-1}x_kln\frac{\mu_k}{1-\sum_{k=1}^{M-1}\mu_k}+ln(1-\sum_{k=1}^{M-1}\mu_k)\}\\ &=(1-\sum_{k=1}^{M-1}\mu_k)exp\{\sum_{k=1}^{M-1}x_kln\frac{\mu_k}{1-\sum_{k=1}^{M-1}\mu_k}\}\\ \end{aligned}
  
  所以 natural parameter 正是：
  
  ηk=lnμk1−∑M−1j=1μj
  
  \eta_k=ln\frac{\mu_k}{1-\sum_{j=1}^{M-1}\mu_j}
  这里的 η\eta 是一个向量，比 Bernoulli 要复杂，因为需要考虑 MM 个不同分类。上述公式中的分母其实就是 M−1M-1 以外的那一个分类的概率 μk=M\mu_{k=M}，所以其实也有点 odds 的意思；这里可以理解为我们随意选择了一个分类作为base，然后用其他分类出现的概率对其求对数比例，把可能性的取值范围扩展到了 (−∞,+∞)(-\infty,+\infty)。 作为被选择作base的分类，其 ηk=M=ln1=0\eta_{k=M}=ln1=0 。
  
  下面推导其 inverse function 即 GLM 使用的 response function，这个过程比 logistic 要复杂很多。首先等价变换 link function：
  
  lnμk1−∑M−1j=1μj=ηkμk1−∑M−1j=1μj=exp{ηk}
  
  ln\frac{\mu_k}{1-\sum_{j=1}^{M-1}\mu_j}=\eta_k\\ \frac{\mu_k}{1-\sum_{j=1}^{M-1}\mu_j}=exp\{\eta_k\}
  接下来，对上式 left side 累加 M 个分类的值：
  
  left=∑k=1Mμk1−∑M−1j=1μj=11−∑M−1j=1μj∑k=1Mμk=11−∑M−1j=1μj
  
  left=\sum^{M}_{k=1}\frac{\mu_k}{1-\sum_{j=1}^{M-1}\mu_j}=\frac{1}{1-\sum_{j=1}^{M-1}\mu_j}\sum^{M}_{k=1}\mu_k=\frac{1}{1-\sum_{j=1}^{M-1}\mu_j}
  对 right side 累加 M 个分类的值：
  
  right=∑k=1Mexp{ηk}=exp{ηk=M}+∑j=1M−1exp{ηj}=1+∑j=1M−1exp{ηj}
  
  right=\sum^{M}_{k=1}exp\{\eta_k\}=exp\{\eta_{k=M}\}+\sum^{M-1}_{j=1}exp\{\eta_j\}=1+\sum^{M-1}_{j=1}exp\{\eta_j\}
  
  俩个式子结合则有：
  
  11−∑M−1j=1μj=1+∑j=1M−1exp{ηj}
  
  \frac{1}{1-\sum_{j=1}^{M-1}\mu_j}=1+\sum^{M-1}_{j=1}exp\{\eta_j\}
  重新代入 link function 则有：
  
  ηk=lnμk1−∑M−1j=1μj=lnμk(1+∑j=1M−1exp{ηj})exp{ηk}=μk(1+∑j=1M−1exp{ηj})μk=exp{ηk}1+∑M−1j=1exp{ηj}
  
  \eta_k=ln\frac{\mu_k}{1-\sum_{j=1}^{M-1}\mu_j}=ln\mu_k(1+\sum^{M-1}_{j=1}exp\{\eta_j\})\\ exp\{\eta_k\}=\mu_k(1+\sum^{M-1}_{j=1}exp\{\eta_j\})\\ \mu_k=\frac{exp\{\eta_k\}}{1+\sum^{M-1}_{j=1}exp\{\eta_j\}}
  这里对于特殊的 μk=M=exp{ηk=M}1+∑M−1j=1exp{ηj}=11+∑M−1j=1exp{ηj}\mu_{k=M}=\displaystyle\frac{exp\{\eta_{k=M}\}}{1+\sum^{M-1}_{j=1}exp\{\eta_j\}}=\displaystyle\frac{1}{1+\sum^{M-1}_{j=1}exp\{\eta_j\}}，所以对任意分类上式都是成立的。
  最终，softmax 的形式为：
  
  μk=exp{ηk}1+∑M−1j=1exp{ηj}
  
  \mu_k=\frac{exp\{\eta_k\}}{1+\sum^{M-1}_{j=1}exp\{\eta_j\}}
  也可以等价地写作：
  
  μk=exp{ηk}∑Mj=1exp{ηj}
  
  \mu_k=\frac{exp\{\eta_k\}}{\sum^{M}_{j=1}exp\{\eta_j\}}
  这个结果近似于 logistic/sigmoid 的：
  
  μ=11+exp{−η}=exp{η}1+exp{η}
  
  \mu=\frac{1}{1+exp\{-\eta\}}=\frac{exp\{\eta\}}{1+exp\{\eta\}}
  且 logistic/sigmoid 中第二种可能的概率为：
  
  1−μ=1−exp{η}1+exp{η}=11+exp{η}
  
  1-\mu=1-\frac{exp\{\eta\}}{1+exp\{\eta\}}=\frac{1}{1+exp\{\eta\}}
  可见 logistic/sigmoid 只是 softmax 在 M=2M=2 时的特殊形式。