20. 多项式模型¶

20.1. 类别分布¶

类别分布可以看做是伯努利分布的扩展，在讲类别分布之前，先回顾一下伯努利分布。伯努利变量只有两个离散状态，一般用 \(0\) 和 \(1\) 表示，伯努利分布的概率质量函数为

(20.1.1)¶\[f(y;\pi) = \pi^y (1-\pi)^{(1-y)}\]

为方便扩展到类别分布，我们先定义一个指示函数， \(\mathbb{I}(true)=1\) ，指示函数只有两个值 \(0 \text{或者} 1\) ，当满足条件时函数值为 \(1\) ，否则为 \(0\) 。

(20.1.2)¶\[\begin{split}\mathbb{I}(x=A)=\left\{ \begin{aligned} 1 & \quad &if \quad x = A \\ 0 & \quad &otherwise \end{aligned} \right.\end{split}\]

用指示函数可以把公式(20.1.1) 改写成

(20.1.3)¶\[f(y;\pi) = \prod_{k} \pi_k^{\mathbb{I}(y=k)}, \quad k \in \{0,1\}\]

其中 \(k\) 表示伯努利分布的两个状态，\(k=0\) 表示状态 \(0\) ，\(k=0\) 表示状态 \(1\) 。 \(\pi_0\) 表示变量 \(Y=0\) 的概率，\(\pi_1\) 表示变量 \(Y=1\) 的概率。。

(20.1.4)¶\[ \begin{align}\begin{aligned}P(Y=0) = f(y=0;\pi) = \pi_0^1 \pi_1^0=\pi_0\\P(Y=1) = f(y=1;\pi) = \pi_0^0 \pi_1^1=\pi_1\end{aligned}\end{align} \]

参数 \(\pi=[\pi_0,\pi_1]\) 表示的是概率值，满足约束 \(\sum_k \pi_k = \pi_0 + \pi_1 =1\) ，\(\pi_0\) 可以由 \(\pi_1\) 计算得到 \(\pi_0=1-\pi_1\) 。因此对于伯努利分布来说，实际上不需要两个未知参数，只需要一个参数就可以，伯努利分布的概率质量函数一般会写成公式(20.1.1) 的单参数形式。

伯努利变量只有两个状态，如果增加到多个状态就是类别分布。类别分布是伯努利分布的扩展，由两个离散状态扩展到更多的离散状态。假设随机变量 \(Y\) 是 \(K\) 个离散状态的类别变量，公式(20.1.1) 可以扩展成类别分布的概率质量函数。

(20.1.5)¶\[ \begin{align}\begin{aligned}f(y;\pi_1,\pi_2,\dots,\pi_K) &= \pi_1^{\mathbb{I}(y=1)} \pi_2^{\mathbb{I}(y=2)} \cdots \pi_K^{\mathbb{I}(y=K)}\\&= \prod_{k=1}^{K} \pi_k^{\mathbb{I}(y=k)} \quad ( 1 \le k \le K)\end{aligned}\end{align} \]

\(\pi_k\) 表示类别 \(k\) 的概率值，\(P(Y=y_k)=\pi_k\) ，所有的 \(\pi_k\) 满足概率约束

(20.1.6)¶\[\sum_{k=1}^K \pi_k = 1\]

既然所有的 \(\pi_k\) 相加为 \(1\) ，就没有必要使用 \(K\) 个参数来表示 \(K\) 个类别的概率，最后一个类别 \(K\) 的概率可以写成

(20.1.7)¶\[\pi_K = 1 - \sum_{k=1}^{K-1} \pi_k\]

公式(20.1.5) 形式有 \(K\) 个参数，可以转换成只有 \(K-1\) 个参数的形式。

(20.1.8)¶\[ \begin{align}\begin{aligned}f(y;\pi_1,\pi_2,\dots,\pi_{K-1}) &= \pi_1^{\mathbb{I}(y=1)} \pi_2^{\mathbb{I}(y=2)} \cdots \pi_{(K-1)}^{\mathbb{I}(y=K-1)} (1 - \sum_{k=1}^{K-1} \pi_k)^{\mathbb{I}(y=K)}\\&= (1 - \sum_{k=1}^{K-1} \pi_k)^{\mathbb{I}(y=K)} \prod_{k=1}^{K-1} \pi_k^{\mathbb{I}(y=k)}\end{aligned}\end{align} \]

由于类别变量 \(Y\) 有多个离散状态值，其期望是一个向量，方差是一个协方差矩阵。类别分布的期望向量为

(20.1.9)¶\[\mathbb{E}[Y] = [\pi_1,\pi_2,\dots,\pi_K]^T\]

方差 \(V(Y_k)\) 和协方差 \(Cov(Y_p,Y_q)\) 为

(20.1.10)¶\[ \begin{align}\begin{aligned}V(Y_k) &= \pi_k(1-\pi_k)\\Cov(Y_p,Y_q) &= - \pi_p \pi_q, \quad (p \ne q)\end{aligned}\end{align} \]

重要

注意，类别变量 \(Y=[y_1,y_2,\dots,y_K]\) 的各个离散状态值是没有任何大小和顺序的关系，各个离散状态值是不可比较的。

20.2. softmax 回归模型¶

20.2.1. 模型定义¶

现在我们推导下 GLM 中类别分布的回归模型，在推导过程中，需要用如下一个等式。

(20.2.1)¶\[{\mathbb{I}(y=K)} = 1 - \sum_{k=1}^{K-1} {\mathbb{I}(y=k)}\]

类别分布的概率质量函数公式(20.1.5) 转换成指数族的形式过程为

(20.2.2)¶\[ \begin{align}\begin{aligned}f(y;\pi_1,\pi_2,\dots,\pi_{K-1}) &= \pi_K^{\mathbb{I}(y=K)} \prod_{k=1}^{K-1} \pi_k^{\mathbb{I}(y=k)}\\&= \exp \left [ \ln \left ( \pi_K^{\mathbb{I}(y=K)} \right )+ \ln \left ( \prod_{k=1}^{K-1} \pi_k^{\mathbb{I}(y=k)} \right ) \right ]\\&= \exp \left [ {\mathbb{I}(y=K)} \ln \pi_K + \sum_{k=1}^{K-1} \ln \pi_k^{{\mathbb{I}(y=k)}} \right ]\\&= \exp \left \{ \left [ 1 - \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} \right ] \ln \pi_K + \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} \ln \pi_k \right \}\\ &= \exp \left \{ \ln \pi_K - \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} \ln \pi_K + \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} \ln \pi_k \right \}\\ &= \exp \left \{ \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} \ln \frac{\pi_k}{\pi_K} +\ln \pi_K \right \}\end{aligned}\end{align} \]

上式中存在指示函数不方便处理，为方便处理需要转换一下。响应变量 \(Y\) 是一个有 \(K\) 个可能取值的离散变量，我们把 \(Y\) 转换成一个长度为 \(K\) 的向量。令 \(T(Y)\) 表示一个长度为 \(K\) 的向量，向量 \(T(Y)\) 只能有一个元素为 \(1\) ，其余元素为 \(0\) ，\(T(Y)_k=1\) 表示向量 \(T(Y)\) 第 \(k\) 个元素为 \(1\) ，等价于响应变量 \(Y=k\) 。

(20.2.3)¶\[\begin{split}T(y)_1 = \begin{bmatrix} 1 \\ 0 \\0 \\ \vdots \\0 \\ 0 \end{bmatrix}, T(y)_2 = \begin{bmatrix} 0 \\ 1 \\0 \\ \vdots \\0 \\ 0 \end{bmatrix}, \dots, T(y)_{K-1} = \begin{bmatrix} 0 \\ 0 \\0 \\ \vdots \\1 \\ 0 \end{bmatrix}, T(y)_{K} = \begin{bmatrix} 0 \\ 0 \\0 \\ \vdots \\0 \\ 1 \end{bmatrix}\end{split}\]

指示函数和向量 \(T(Y)\) 的关系可以写为

(20.2.4)¶\[T(Y)_k = \mathbb{I}(y=k)\]

此外，我们再定义一个长度为 \(K\) 的参数向量 \(\theta\) ，

(20.2.5)¶\[\begin{split}\theta = \begin{bmatrix} \ln \pi_1/\pi_K \\ \ln (\pi_2/\pi_K) \\ \vdots \\ \ln (\pi_k/\pi_K) \\ \vdots \\ \ln (\pi_{K-1}/\pi_K) \\ \ln (\pi_{K}/\pi_K)=0 \end{bmatrix}\end{split}\]

依据向量內积的定理则有

(20.2.6)¶\[ \begin{align}\begin{aligned}\theta^T T(y) &= \sum_{k=1}^{K} T(y)_k \ln \frac{\pi_k}{\pi_K}\\&= T(y)_K \ln \frac{\pi_K}{\pi_K} + \sum_{k=1}^{K-1} T(y)_k \ln \frac{\pi_k}{\pi_K}\\&= 0 + \sum_{k=1}^{K-1} T(y)_k \ln \frac{\pi_k}{\pi_K}\\&= \sum_{k=1}^{K-1} T(y)_k \ln \frac{\pi_k}{\pi_K}\\&= \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} \ln \frac{\pi_k}{\pi_K}\end{aligned}\end{align} \]

此外根据概率约束有

(20.2.7)¶\[\pi_K = 1 - \sum_{k=1}^{K-1} \pi_k\]

通过把累加转化成向量內积的方式，把公式(20.2.2) 写成指数族的形式

(20.2.8)¶\[ \begin{align}\begin{aligned}f(y;\theta) &= \exp \left \{ \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} \ln \frac{\pi_k}{ \pi_K} + \ln \pi_K \right \}\\&= \exp \left \{ \theta^T T(y) + A(\theta) \right \}\end{aligned}\end{align} \]

上式就是类别分布概率质量函数的指数族形式，其中 \(\theta\) 就是自然参数，并且 \(\theta\) 和 \(T(y)\) 分别是一个向量。可以看到类别分布的指数族形式是不存在分散函数 \(a(\phi)\) 的，或者说 \(a(\phi)=1\) 。

标准连接函数

自然参数 \(\theta\) 和期望参数 \(\mu=\pi\) 的关系为

(20.2.9)¶\[\theta_k = \ln \frac{ \pi_k}{ \pi_K} = \ln \frac{ \pi_k}{ ( 1 - \sum_{k=1}^{K-1} \pi_k) }\]

根据标准连接函数的定义，当线性预测 \(\eta\) 等于自然参数 \(\theta\) 时得到的就是标准连接函数因此其标准连接函数为

(20.2.10)¶\[ \begin{align}\begin{aligned}\eta_k = \theta_k &= \ln \frac{\pi_k}{\pi_K}\\&= \ln \frac{ \pi_k}{ ( 1 - \sum_{k=1}^{K-1} \pi_k) }\\&= \ln \frac{ \mu_k}{ ( 1 - \sum_{k=1}^{K-1} \mu_k) }\end{aligned}\end{align} \]

响应函数

响应函数是连接函数的反函数，现在我们推导下公式(20.2.10) 的反函数。首先等号两边取值指数操作，可得

(20.2.11)¶\[\pi_k = \pi_K e^{\eta_k}\]

上式中仍然存在 \(\pi_K\) ，需要推导下 \(\pi_K\) 如何用 \(\eta\) 表示。

(20.2.12)¶\[ \begin{align}\begin{aligned}1 &= \sum_{k=1}^K \pi_k\\&= \pi_K + \sum_{k=1}^{K-1} \pi_k\\&= \pi_K + \sum_{k=1}^{K-1} \pi_K e^{\eta_k}\\&= \pi_K (1+ \sum_{k=1}^{K-1} e^{\eta_k})\end{aligned}\end{align} \]

因此，有

(20.2.13)¶\[\pi_K = \frac{1}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}}\]

代入到公式(20.2.10) 可得

(20.2.14)¶\[\mu_k = \pi_k = \pi_K e^{\eta_k} = \frac{e^{\eta_k}}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}} , \quad for \quad k=1,2,\dots,K-1\]

公式(20.2.14) 就是类别分布的响应函数，此函数通常叫作 softmax 函数，它是 logistic 函数的扩展，当 \(K=2\) 时，softmax 就是退化成 logistic 函数。

对于有 \(K\) 个类别的 softmax 回归模型来说，其线性预测器只有 \(K-1\) 个，这意味着有 \(K-1\) 个协变量参数向量， \(\beta=[\beta_1,\beta_2,\dots,\beta_{K-1}]\) 。

(20.2.15)¶\[\eta_k = \ln \frac{ \mu_k}{ ( 1 - \sum_{k=1}^{K-1} \mu_k) } = \beta_k x^T , \quad for \quad k=1,2,\dots,K-1\]

这 \(K-1\) 个类别的响应函数就是 softmax 函数。

(20.2.16)¶\[\mu_k = \frac{e^{\eta_k}}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}} , \quad for \quad k=1,2,\dots,K-1\]

第 \(K\) 个类别不需要独立的线性预测器，依据概率的约束，它的期望可以通过其余 \(K-1\) 个计算得到。

(20.2.17)¶\[\mu_K = 1- \sum_{k=1}^{K-1} \mu_k = \frac{1}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}}\]

最后我们整理下 softmax 回归的各个关键部分。

(20.2.18)¶\[ \begin{align}\begin{aligned}\text{自然参数} & \quad \theta_k = \ln \frac{ \pi_k}{ \pi_K} = \ln \frac{ \pi_k}{ ( 1 - \sum_{k=1}^{K-1} \pi_k) }\\ \text{累积分布函数} & \quad b(\theta) = \ln \pi_K = ( 1 - \sum_{k=1}^{K-1} \pi_k)\\\text{期望} & \quad \mu_k = b'(\theta_k) = \pi_k\\\text{方差} & \quad V(Y_k) = \pi_k(1-\pi_k)\\\text{协方差} & \quad Cov(Y_p,Y_q) = - \pi_p \pi_q, \quad (p \ne q)\\\text{分散函数} & \quad a(\phi) = 1\\\text{标准连接函数} & \quad g(\mu_k) = \ln \frac{ \pi_k}{ \pi_K } = \ln \frac{ \pi_k}{ ( 1 - \sum_{k=1}^{K-1} \pi_k) }\\\text{标准连接函数一阶导} & \quad g'(\mu_k) = \frac{\pi_k + \pi_K}{\pi_k \pi_K} =\frac{\pi_k+1-\sum_{k=1}^{K-1} \pi_k}{\pi_k( 1 - \sum_{k=1}^{K-1} \pi_k)}\\\text{响应函数} & \quad r(\eta_k) = \frac{e^{\eta_k}}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}}\end{aligned}\end{align} \]

类别分布是伯努利分布在多类别上的扩展，因此二者在连接函数、响应函数等方面都有一定的关联性。二值离散变量的概率分布是伯努利分布，其标准连接对应的响应函数是 logistic 函数，其回归模型是称为 logistic 回归，可用于处理二分类响应数据。多值离散变量的概率分布是类别分布，其标准连接对应的响应函数是 softmax 函数，其回归模型称为 softmax 回归，可用于处理多分类的响应数据。由于 softmax 是 logistic 的扩展，有时也会把 softmax 回归称为多元逻辑回归。

20.2.2. 参数估计¶

softmax 回归模型也是 GLM 中的一员，并且它是 logistic 回归模型的扩展，它在 IRLS 估计算法的各个组件和 logistic 模型是类似的。

softmax 回归模型的对数似然函数为

(20.2.19)¶\[ \begin{align}\begin{aligned}\ell(\hat{\pi};y) &= \sum_{i=1}^N \left \{ \sum_{k=1}^{K-1} \left [ \mathbb{I}(y_i=k) \ln \frac{ \hat{\pi}_{ik}}{\pi_{iK}} \right ] + \ln \hat{\pi}_{iK} \right \}\\&= \sum_{i=1}^N \left \{ \sum_{k=1}^{K-1} \left [ T(y_i)_k \ln \frac{\pi_{ik}}{\hat{\pi}_{iK}} \right ] + \ln \hat{\pi}_{iK} \right \}\\&= \sum_{i=1}^N \left \{ \sum_{k=1}^{K-1} \left [ T(y_i)_k \ln \hat{\pi}_{ik} \right ] - \sum_{k=1}^{K-1} \left [ T(y_i)_k \ln {\hat{\pi}_{iK}} \right ] + \ln \hat{\pi}_{iK} \right \}\\ &= \sum_{i=1}^N \left \{ \sum_{k=1}^{K-1} \left [ T(y_i)_k \ln \hat{\pi}_{ik} \right ] - \sum_{k=1}^{K-1} \left [ T(y_i)_k \ln {\hat{\pi}_{iK}} \right ] + \ln \hat{\pi}_{iK} \right \}\end{aligned}\end{align} \]

不同于二元逻辑回归模型，softmax 回归模型有 \(K-1\) 个协变量参数向量，参数向量 \(\beta_k\) 的权重矩阵 \(W_k\) 和工作响应矩阵 \(Z_k\) 的计算公式分别为：

(20.2.20)¶\[ \begin{align}\begin{aligned}W_{iik} &= \frac{ 1}{ a(\phi) \nu(\hat{\mu}_{ik}) ( g_{ik}' )^2}\\ &= \frac{ \hat{\pi}_{ik} \hat{\pi}_{iK}^2 }{ (\hat{\pi}_{ik} + \hat{\pi}_{iK})^2(1-\hat{\pi}_{ik}) }\end{aligned}\end{align} \]

(20.2.21)¶\[ \begin{align}\begin{aligned}Z_{ik} &= [T(y_i)_{k}- \hat{\mu}_{ik}] g_{ik}' + \eta_{ik}\\ &= \frac{[T(y_i)_{k}- \hat{\mu}_{ik}](\hat{\pi}_{ik} +\hat{\pi}_{iK}) } {\hat{\pi}_{ik} \hat{\pi}_{iK}} + \eta_{ik}\end{aligned}\end{align} \]

20.3. 多项式分布¶

伯努利分布表示单次伯努利实验成功或者失败的概率，多次伯努利实验成功次数的概率分布是二项式分布。按照这个方式，多次类别分布实验各个状态发生次数的概率分布就叫做多项式分布(multinomial distribution)。假设实验总次数为 \(n\) ，每个类别发生的次数为 \(y_k\) ，多项式分布的概率质量函数为

(20.3.1)¶\[f(y;\pi) = \frac{n!}{y_1!y_2!\dots y_K! }\prod_{k=1}^{K} \pi_k^{y_k}\]

多项式分布变量的期望为

(20.3.2)¶\[\mathbb{E}[y_k] = n \pi_k\]

方差 \(V(Y_k)\) 和协方差 \(Cov(Y_p,Y_q)\) 为

(20.3.3)¶\[ \begin{align}\begin{aligned}V(Y_k) &= n \pi_k(1-\pi_k)\\Cov(Y_p,Y_q) &= - n \pi_p \pi_q, \quad (p \ne q)\end{aligned}\end{align} \]

可以看到多项式分布和类别分布的区别就是多了一个实验的总次数 \(n\) ，当 \(n=1\) 时，多项式分布就等价于类别分布，类别分布是多项式分布的一个特例。当 \(K=2\) 时，多项式分布就等价于二项式分布，因此二项式分布也是多项式分布的一个特例。多项式分布于类别分布的关系，就好比二项式分布于伯努利分布的关系，读者可以回顾一下二项式模型的章节加深理解。此外多项式分布也是二项式分布的扩展，二项式分布变量只有两个状态，而多项式分布扩展的更多的状态。

需要注意的是，虽然通常会使用连续的整数表示多项式变量的各个离散状态，但这些状态间是没有大小顺序关系的。存在顺序关系的离散状态随机变量，我们下一章再探讨。

20.4. 多项式回归模型¶

类别分布可以用来处理多分类的响应数据，而多项式分布分科

(20.4.1)¶\[ \begin{align}\begin{aligned}f(y;\pi) &= \frac{n!}{y_1!y_2!\dots y_K! }\prod_{k=1}^{K} \pi_k^{y_k}\\&=\exp \left \{ \ln \left [ \frac{n!}{y_1!y_2!\dots y_K! } \right ] + \ln \prod_{k=1}^{K} \pi_k^{y_k} \right \}\\&= \exp \left \{ \sum_{k=1}^{K} {y_k} \ln \pi_k + \ln \left [ \frac{n!}{y_1!y_2!\dots y_K! } \right ] \right \}\end{aligned}\end{align} \]