##################################### 多项式模型 ##################################### https://www.freesion.com/article/4717201261/ https://www.freesion.com/article/4717201261/ https://www.jianshu.com/p/99c4238f067e 类别分布 ##################################### 类别分布可以看做是伯努利分布的扩展，在讲类别分布之前，先回顾一下伯努利分布。 伯努利变量只有两个离散状态，一般用 :math:`0` 和 :math:`1` 表示，伯努利分布的概率质量函数为 .. math:: :label: eq_mg_003 f(y;\pi) = \pi^y (1-\pi)^{(1-y)} 为方便扩展到类别分布，我们先定义一个指示函数， :math:`\mathbb{I}(true)=1` ，指示函数只有两个值 :math:`0 \text{或者} 1` ，当满足条件时函数值为 :math:`1` ，否则为 :math:`0` 。 .. math:: \mathbb{I}(x=A)=\left\{ \begin{aligned} 1 & \quad &if \quad x = A \\ 0 & \quad &otherwise \end{aligned} \right. 用指示函数可以把 :eq:`eq_mg_003` 改写成 .. math:: f(y;\pi) = \prod_{k} \pi_k^{\mathbb{I}(y=k)}, \quad k \in \{0,1\} 其中 :math:`k` 表示伯努利分布的两个状态，:math:`k=0` 表示状态 :math:`0` ，:math:`k=0` 表示状态 :math:`1` 。 :math:`\pi_0` 表示变量 :math:`Y=0` 的概率 ，:math:`\pi_1` 表示变量 :math:`Y=1` 的概率。 。 .. math:: 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 参数 :math:`\pi=[\pi_0,\pi_1]` 表示的是概率值，满足约束 :math:`\sum_k \pi_k = \pi_0 + \pi_1 =1` ，:math:`\pi_0` 可以由 :math:`\pi_1` 计算得到 :math:`\pi_0=1-\pi_1` 。因此对于伯努利分布来说，实际上不需要两个未知参数，只需要一个参数就可以， 伯努利分布的概率质量函数一般会写成 :eq:`eq_mg_003` 的单参数形式。 伯努利变量只有两个状态，如果增加到多个状态就是类别分布。 类别分布是伯努利分布的扩展，由两个离散状态扩展到更多的离散状态。 假设随机变量 :math:`Y` 是 :math:`K` 个离散状态的类别变量， :eq:`eq_mg_003` 可以扩展成类别分布的概率质量函数。 .. math:: :label: eq_multinomial_011 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) :math:`\pi_k` 表示类别 :math:`k` 的概率值，:math:`P(Y=y_k)=\pi_k` ， 所有的 :math:`\pi_k` 满足概率约束 .. math:: :label: eq_multinomial_012 \sum_{k=1}^K \pi_k = 1 既然所有的 :math:`\pi_k` 相加为 :math:`1` ，就没有必要使用 :math:`K` 个参数来表示 :math:`K` 个类别的概率 ，最后一个类别 :math:`K` 的概率可以写成 .. math:: \pi_K = 1 - \sum_{k=1}^{K-1} \pi_k :eq:`eq_multinomial_011` 形式有 :math:`K` 个参数，可以转换成只有 :math:`K-1` 个参数的形式。 .. math:: :label: eq_multinomial_013 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)} 由于类别变量 :math:`Y` 有多个离散状态值，其期望是一个向量，方差是一个协方差矩阵。 类别分布的期望向量为 .. math:: \mathbb{E}[Y] = [\pi_1,\pi_2,\dots,\pi_K]^T 方差 :math:`V(Y_k)` 和协方差 :math:`Cov(Y_p,Y_q)` 为 .. math:: V(Y_k) &= \pi_k(1-\pi_k) Cov(Y_p,Y_q) &= - \pi_p \pi_q, \quad (p \ne q) .. important:: 注意，类别变量 :math:`Y=[y_1,y_2,\dots,y_K]` 的各个离散状态值是没有任何大小和顺序的关系， 各个离散状态值是不可比较的。 softmax 回归模型 ######################################## 模型定义 ================================= 现在我们推导下 ``GLM`` 中类别分布的回归模型， 在推导过程中，需要用如下一个等式。 .. math:: {\mathbb{I}(y=K)} = 1 - \sum_{k=1}^{K-1} {\mathbb{I}(y=k)} 类别分布的概率质量函数 :eq:`eq_multinomial_011` 转换成指数族的形式过程为 .. math:: :label: eq_multinomial_020 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 \} 上式中存在指示函数不方便处理，为方便处理需要转换一下。 响应变量 :math:`Y` 是一个有 :math:`K` 个可能取值的离散变量， 我们把 :math:`Y` 转换成一个长度为 :math:`K` 的向量。 令 :math:`T(Y)` 表示一个长度为 :math:`K` 的向量， 向量 :math:`T(Y)` 只能有一个元素为 :math:`1` ，其余元素为 :math:`0` ，:math:`T(Y)_k=1` 表示向量 :math:`T(Y)` 第 :math:`k` 个元素为 :math:`1` ，等价于响应变量 :math:`Y=k` 。 .. math:: 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} 指示函数和向量 :math:`T(Y)` 的关系可以写为 .. math:: T(Y)_k = \mathbb{I}(y=k) 此外，我们再定义一个长度为 :math:`K` 的参数向量 :math:`\theta` ， .. math:: \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} 依据向量內积的定理则有 .. math:: \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} 此外根据概率约束有 .. math:: \pi_K = 1 - \sum_{k=1}^{K-1} \pi_k 通过把累加转化成向量內积的方式，把 :eq:`eq_multinomial_020` 写成指数族的形式 .. math:: 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 \} 上式就是类别分布概率质量函数的指数族形式， 其中 :math:`\theta` 就是自然参数， 并且 :math:`\theta` 和 :math:`T(y)` 分别是一个向量。 可以看到类别分布的指数族形式是不存在分散函数 :math:`a(\phi)` 的， 或者说 :math:`a(\phi)=1` 。 **标准连接函数** 自然参数 :math:`\theta` 和 期望参数 :math:`\mu=\pi` 的关系为 .. math:: :label: eq_mg_025 \theta_k = \ln \frac{ \pi_k}{ \pi_K} = \ln \frac{ \pi_k}{ ( 1 - \sum_{k=1}^{K-1} \pi_k) } 根据标准连接函数的定义，当线性预测 :math:`\eta` 等于自然参数 :math:`\theta` 时得到的就是标准连接函数 因此其标准连接函数为 .. math:: :label: eq_mg_030 \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) } **响应函数** 响应函数是连接函数的反函数，现在我们推导下 :eq:`eq_mg_030` 的反函数。 首先等号两边取值指数操作，可得 .. math:: \pi_k = \pi_K e^{\eta_k} 上式中仍然存在 :math:`\pi_K` ，需要推导下 :math:`\pi_K` 如何用 :math:`\eta` 表示。 .. math:: 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}) 因此，有 .. math:: \pi_K = \frac{1}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}} 代入到 :eq:`eq_mg_030` 可得 .. math:: :label: eq_mg_035 \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 :eq:`eq_mg_035` 就是类别分布的响应函数，此函数通常叫作 ``softmax`` 函数， 它是 ``logistic`` 函数的扩展， 当 :math:`K=2` 时，``softmax`` 就是退化成 ``logistic`` 函数。 对于有 :math:`K` 个类别的 ``softmax`` 回归模型来说，其线性预测器只有 :math:`K-1` 个，这意味着有 :math:`K-1` 个协变量参数向量， :math:`\beta=[\beta_1,\beta_2,\dots,\beta_{K-1}]` 。 .. math:: \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 这 :math:`K-1` 个类别的响应函数就是 ``softmax`` 函数。 .. math:: :label: eq_mg_037 \mu_k = \frac{e^{\eta_k}}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}} , \quad for \quad k=1,2,\dots,K-1 第 :math:`K` 个类别不需要独立的线性预测器， 依据概率的约束，它的期望可以通过其余 :math:`K-1` 个计算得到。 .. math:: \mu_K = 1- \sum_{k=1}^{K-1} \mu_k = \frac{1}{ 1+ \sum_{k=1}^{K-1} e^{\eta_k}} 最后我们整理下 ``softmax`` 回归的各个关键部分。 .. math:: \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}} 类别分布是伯努利分布在多类别上的扩展，因此二者在连接函数、响应函数等方面都有一定的关联性。 二值离散变量的概率分布是伯努利分布，其标准连接对应的响应函数是 ``logistic`` 函数，其回归模型是称为 ``logistic`` 回归， 可用于处理二分类响应数据。 多值离散变量的概率分布是类别分布， 其标准连接对应的响应函数是 ``softmax`` 函数，其回归模型称为 ``softmax`` 回归， 可用于处理多分类的响应数据。 由于 ``softmax`` 是 ``logistic`` 的扩展，有时也会把 ``softmax`` 回归称为多元逻辑回归。 参数估计 ================================= ``softmax`` 回归模型也是 ``GLM`` 中的一员， 并且它是 ``logistic`` 回归模型的扩展， 它在 ``IRLS`` 估计算法的各个组件和 ``logistic`` 模型是类似的。 ``softmax`` 回归模型的对数似然函数为 .. math:: \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 \} 不同于二元逻辑回归模型，``softmax`` 回归模型有 :math:`K-1` 个协变量参数向量， 参数向量 :math:`\beta_k` 的权重矩阵 :math:`W_k` 和工作响应矩阵 :math:`Z_k` 的计算公式分别为： .. math:: 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}) } .. math:: 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} 多项式分布 ######################################## 伯努利分布表示单次伯努利实验成功或者失败的概率， 多次伯努利实验成功次数的概率分布是二项式分布。 按照这个方式， 多次类别分布实验各个状态发生次数的概率分布就叫做多项式分布(multinomial distribution)。 假设实验总次数为 :math:`n` ，每个类别发生的次数为 :math:`y_k` ，多项式分布的概率质量函数为 .. math:: f(y;\pi) = \frac{n!}{y_1!y_2!\dots y_K! }\prod_{k=1}^{K} \pi_k^{y_k} 多项式分布变量的期望为 .. math:: \mathbb{E}[y_k] = n \pi_k 方差 :math:`V(Y_k)` 和协方差 :math:`Cov(Y_p,Y_q)` 为 .. math:: V(Y_k) &= n \pi_k(1-\pi_k) Cov(Y_p,Y_q) &= - n \pi_p \pi_q, \quad (p \ne q) 可以看到多项式分布和类别分布的区别就是多了一个实验的总次数 :math:`n` ， 当 :math:`n=1` 时，多项式分布就等价于类别分布，类别分布是多项式分布的一个特例。 当 :math:`K=2` 时，多项式分布就等价于二项式分布， 因此二项式分布也是多项式分布的一个特例。 多项式分布于类别分布的关系，就好比二项式分布于伯努利分布的关系，读者可以回顾一下二项式模型的章节加深理解。 此外多项式分布也是二项式分布的扩展，二项式分布变量只有两个状态，而多项式分布扩展的更多的状态。 需要注意的是，虽然通常会使用连续的整数表示多项式变量的各个离散状态， 但这些状态间是没有大小顺序关系的。 存在顺序关系的离散状态随机变量，我们下一章再探讨。 多项式回归模型 ##################################### 类别分布可以用来处理多分类的响应数据，而多项式分布分科 .. math:: 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 \}