###############################################
Gamma 模型
###############################################


指数分布是预测 **下一次** 时间发生等待的时间，
更进一步，如果需要预测第k次事件发生需要等待的时间呢？
这就是 Gamma 分布。
在介绍 Gamma 分布之前，先简单的介绍一下 Gamma 函数。

Gamma 函数
#############################


https://www.probabilitycourse.com/chapter4/4_2_4_Gamma_distribution.php


Gamma 分布
#############################

在之前的章节我们已经讨论了如何从泊松过程推导出指数分布，
Gamma 分布的推导的过程是类似的。
不同的地方在于指数分布是等待下一次事件发生，
而 Gamma 分布是等待第k次事件的发生。

我们用符号 :math:`k` 表示事件发生的次数，
符号 :math:`T` 表示直到 :math:`k` 次事件发生等待的时间，也就是目标随机变量。
:math:`\lambda` 表示泊松过程中事件发生的频率（单位时间内发生的次数）。
:math:`p(T>t;k)` 表示直到第 :math:`k` 次事件发生等待的时间 :math:`T` 大于 :math:`t` 的概率。
用符号 :math:`p(k;T=t)` 表示在泊松过程中 :math:`t` 个时间单元内事件发生 :math:`k` 次的概率。


泊松分布表示单位时间窗口内事件发生次数的概率分布，
泊松分布的概率分布函数为：

.. math::


    p(k) = \frac{\lambda^k}{k!}  e^{-\lambda}


其表示在一个单位时间内事件发生 :math:`k` 次的概率，唯一的参数 :math:`\lambda` 表示单位时间内事件发生次数的平均值，
也就是 :math:`k` 的期望值，:math:`\mathbb{E}[x]=\lambda` 。
现在要想计算 :math:`t` 个单位时间内事件发生的次数，只需要用 :math:`\lambda t` 替换上式中的 :math:`\lambda` 即可，
:math:`t` 个时间单元内事件发生 :math:`k` 次的概率为：



.. math::

    p(k;t) = \frac{(\lambda t)^k}{k!}  e^{-\lambda t}


.. note::
    为什么用 :math:`\lambda t` 替换 :math:`\lambda` 就可以？

    泊松分布表示单位时间内，随机事件发生 :math:`k` 次的概率分布。
    这个单位时间并没有具体的长度限制，只要保证每个时间片段的长度相同就可以。
    泊松分布唯一的参数就是每个时间片段内随机事件发生次数的平均值 :math:`\lambda` ，
    如果要算 :math:`t` 个时间片段内随机事件发生次数的概率分布，相当于原来的时间片段扩大了 :math:`t` 倍，
    这个 :math:`t` 倍的时间片段也可以看做是一个"单位时间"，这个新的"单位时间"组成一个新的泊松分布，
    并且其平均值参数 :math:`\lambda_{new}` 也是原来的 :math:`t` 倍，即 :math:`\lambda_{new} =  \lambda t`。

    .. math::

        p(k) = \frac{(\lambda_{new})^k}{k!}  e^{-\lambda_{new}} = \frac{(\lambda t)^k}{k!}  e^{-\lambda t}


泊松分布表示的是 *一个单位时间片段内* 随机事件发生 **次数** 的概率分布，
指数分布表示的是随机事件 *第一次* 发生的需要的 **时间** 的概率分布。
泊松分布描述的是 **次数**，是离散变量的分布；
指数分布描述的是 **时间** ，是连续值变量的分布。
然而指数分布仅仅描述了事件首次发生，不能表示更多次事件发生，不具备一般化，
Gamma 分布就是指数分布的扩展，能够表达事件发生任意次数所需 **时间** 的概率分布。


我们令符号 :math:`p(T>t;k)` 表示随机事件发生 :math:`k` 所需时间大于 :math:`t` 的概率。
如果随机事件第 :math:`k` 发生的时间大于 :math:`t` ，
意味着在 :math:`t` 个时间单元内，事件发生次数一定是小于等于 :math:`k-1` 次。
换句话说，:math:`t` 个时间单元内最多只发生 :math:`k-1` 次，
第 :math:`k` 次发生的时间一定是大于 :math:`t` 。
因此 :math:`p(T>t;k)` 可以看成是在 :math:`t` 个时间单元内，事件发生 :math:`0,1,2,...,k-1` 次的概率之和。


.. math::

    p(T>t;k) = \sum_{i=0}^{k-1} p(i;t)
    = \sum_{i=0}^{k-1} \frac{(\lambda t)^{i}}{i!}  e^{-\lambda t}


现在我们得到了第 :math:`k` 次事件发生时间大于 :math:`t` 的概率 :math:`p(T>t;k)`
，反过来，第 :math:`k` 次事件发生时间小于等于 :math:`t` 的概率为：


.. math::
    :label: eq_gamma_008

    p(T\le t;k) &= 1-p(T>t;k)

    &= 1- \sum_{i=0}^{k-1} \frac{(\lambda t)^{i}}{i!}  e^{-\lambda t}

显然 :eq:`eq_gamma_008` 是一个累积分布函数(Cumulative Distribution Function)，
是概率密度函数的积分，对其进行微分可以得到概率密度函数，
这里省略微分的过程，直接给出结果


.. math::
    :label: eq_gamma_009

    f(t)
    &= \frac{\lambda e^{-\lambda t}(\lambda t)^{k-1}}{(k-1)!}

    &= \frac{\lambda^k e^{-\lambda t}t^{k-1}}{\Gamma(k)}

:eq:`eq_gamma_009` 就是 Gamma 分布的概率密度函数，其表示随机事件发生 :math:`k` 次所需时间 :math:`t`
的概率分布。注意，:math:`f(t)` 不是具体的时间值，而是时间 :math:`t` 的概率分布。
参数 :math:`k` 表示事件发生的次数，又被称为形状参数(shape parameter)。
参数 :math:`\lambda` 来源于泊松分布，表示随机事件发生的频次，即单位时间发生的平均次数，又被称为速率参数(rate parameter)。


.. note::
    什么是形状参数(shape parameter)？

    在概率分布中，按照参数对概率分布函数曲线的影响，可以分为几种。

    - 位置参数(location parameter)，影响着图形在 :math:`x` 轴上的位置；
    - 尺度参数(scale parameter)，控制着图形的拉伸和缩小，可以缩放图形。
      尺度参数的倒数称为速率参数(rate parameter)，对图形的影响与尺度参数是一样。
    - 形状参数(shape parameter)，其既不是位置参数也不是尺度参数（也不是关于这两者的函数）。
      形状参数直接影响图形分布的形状，而不是简单地移动分布（如位置参数）或拉伸/缩小分布（如比例参数）。


Gamma 分布的期望和方差分别为：

.. math::

    \mathbb{E}[T] &= \frac{k}{\lambda}

    Var(T) &=  \frac{k}{\lambda^2}

:math:`k` 是事件发生总次数，:math:`\lambda` 是事件发生速率，
显然事件发生 :math:`k` 次所需要的平均时间就是 :math:`k/\lambda` ，这和 Gamma 分布的期望是一致的。

当 :math:`k=1` 时，Gamma分布就退化为指数分布:

.. math::

    f(t;k=1,\lambda) =  \lambda e^{-\lambda t} \quad  (\lambda ,t > 0)


因此有 :math:`Gamma(1,\lambda) = Exponential(\lambda)`
成立，更一般的，如果有 :math:`n` 个独立的指数分布 :math:`Exponential(\lambda)` 随机变量，
就可以得到一个 Gamma 分布的随机变量 :math:`Gamma(n,\lambda)`
。

现在我们来看下形状参数和速率参数分别对概率分布函数的影响是怎样的。
:numref:`fg_gamma_001` 展示了形状参数 :math:`k` 对概率分布函数的影响，
:numref:`fg_gamma_002` 展示了速率参数 :math:`\lambda` 对概率分布函数的影响。



.. _fg_gamma_001:

.. figure:: pictures/gamma_k.png
    :scale: 70 %
    :align: center

    Gamma分布不同 :math:`k` 值下图形比较

.. _fg_gamma_002:

.. figure:: pictures/gamma_lambda.png
    :scale: 70 %
    :align: center

    Gamma分布不同 :math:`\lambda` 值下图形比较


通常 Gamma 分布的概率密度函有多种参数化方式，
在计量经济学和其它一些自然科学领域，经常使用形状参数和尺度参数进行参数化表示。

**常见参数化方式1：**

令参数 :math:`\alpha=k` 表式形状参数，参数 :math:`\beta=1/\lambda` 表示尺度参数，
服从 Gamma 分布的随机变量 :math:`X` 的概率密度函数为

.. math::
    :label: eq_gamma_011

    f(x;\alpha,\beta) = \frac{ x^{\alpha-1} e^{-\frac{x}{\beta} }}{\beta^{\alpha} \Gamma(\alpha)}
    \quad (x,\alpha,\beta > 0)

这种形式常见于计量经济学和其它一些自然科学领域。
这种形式下，``Gamma`` 分布的期望和方差为别为

.. math::

    \mathbb{E}[X] &= \alpha \beta

    Var(X) &=  \alpha \beta^2

**常见参数化方式2：**

令参数 :math:`\alpha=k` 表式形状参数，参数 :math:`\beta=\lambda` 表示尺度参数，
这种方式下， 尺度参数 :math:`\beta` 和 :eq:`eq_gamma_011` 是倒数关系，这只是不同的参数化方法而已，
二者是等价的，这时 ``Gamma`` 分布的概率分布函数写成

.. math::
    :label: eq_gamma_012

    f(x;\alpha,\beta) = \frac{ x^{\alpha-1} e^{-\beta x  }  \beta^{\alpha} }{ \Gamma(\alpha)}
    \quad (x,\alpha,\beta > 0)

此时的期望和方差为别为

.. math::

    \mathbb{E}[X] &= \frac{\alpha}{\beta}

    Var(X) &=  \frac{\alpha} {\beta^2}



Gamma 回归模型
####################################




.. n理想情况下，Gamma模型最好与具有恒定变化系数的正响应一起使用。
    但是，该模型对于与后一个准则的较大偏差具有鲁棒性。
    因为两参数Gamma分布的形状是灵活的，
    并且可以参数化以适合许多响应数据形状，
    所以对于许多严格的正响应数据情况，
    它可能比高斯模型更好用。


.. x传统的GLM框架限制模型只能有一个参数：均值 :math:`\mu` ，
    当方差或者连接函数中包含分散参数 :math:`\phi` 时，通常会模型限制其为常量1，或者由用户指定一个值。
    但是，在稍后的内容中，不再是直接的GLM模型而是拟似然(quasi-likelihood)模型。
    后一种情况的主要示例是负二项式模型的传统GLM版本， 前一种情况的示例包括高斯模型，伽马模型和高斯逆模型。



.. x尽管传统的迭代加权最小二乘（IRLS）GLM算法没有像FIML算法那样正式估计分散参数或辅助参数，
    但是分散参数可以从基于IRLS Pearson的离散度或分散统计量估算。
    也就是说，在使用FIML的伽马回归模型中，分散参数的值与基于IRLS Pearson的估算值几乎相同。


.. x在本章中，我们将讨论gamma模型的各种GLM的参数化方法。
    但是，我们将不会讨论三参数模型，也不会对包含审查的伽马模型进行实质性讨论。
    GLM Gamma模型可用于直接对指数回归建模，但我们将这些讨论类型放到本章的最后部分。

在GLM中，Gamma模型用于响应数据只能取大于或等于0的连续值数据进行建模。


Gamma 分布是是一个双参数分布，包含形状参数 :math:`\alpha` 和尺度化参数 :math:`\beta`
，形状参数 :math:`\alpha` 表示事件的发生次数，尺度化参数 :math:`\beta` 表示平均一次事件发生需要的时间。
二者的乘积 :math:`\mu=\alpha \beta` 就是事件发生  :math:`\alpha` 次所需的平均时间，
即 Gamma 分布的期望值(均值)。

当把 Gamma 分布作为 GLM 家族的成员时，
需要把 :math:`\alpha` 看做一个已知的常量，
也就是人为给 :math:`\alpha` 设置一个常数，并且对于所有观测样本都是一样的值。
通常这个值一个经验值，需要根据观测数据分布情况


GLM中指数族分布的标准形式为

.. math::
    :label: eq_gamma_012_1

    f(y;\theta,\phi) = \exp \left \{\frac{\theta y - b(\theta)}{a(\phi)} + c(y,\phi) \right \}


现在我们把 Gamma 分布的概率密度函数转化成上述指数族分布的标准形式，
首先，从上文已知 Gamma 分布的期望为 :math:`\mu=\alpha \beta`
。令 :math:`\phi = 1/\alpha` ，则有 :math:`\alpha = 1 / \phi ,\beta=\mu / \alpha=\mu \phi`
，代入到 :eq:`eq_gamma_011` 可得



.. math::
    :label: eq_gamma_013

    f(y;\mu,\phi) &=  \frac{ y^{\alpha-1} e^{-\frac{y}{\beta} }}{\beta^{\alpha} \Gamma(\alpha)}

    &= \exp \left [
        -\frac{y}{\beta} + (\alpha-1) \ln y - \alpha \ln \beta  - \ln \Gamma(\alpha)
    \right ]

    &= \exp \left [
        -y / \beta  - \alpha \ln \beta  + (\alpha-1) \ln y   - \ln \Gamma(\alpha)
    \right ]


    &= \exp \left [
        \frac{-y}{\mu \phi}  - \frac{ \ln (\mu \phi) }{ \phi} + (\frac{1}{\phi}-1) \ln y   - \ln [\Gamma(1 / \phi)]
    \right ]

    &= \exp \left [
        \frac{y(1/\mu)}{-\phi}  - \frac{ \ln \mu }{ \phi} - \frac{ \ln  \phi }{ \phi}  + (\frac{1-\phi}{\phi}) \ln y  - \ln [\Gamma(1 / \phi)]
    \right ]

    &= \exp \left [
        \frac{y(1/\mu)}{-\phi}  + \frac{ \ln \mu }{ -\phi} - \frac{ \ln  \phi }{ \phi}  + (\frac{1-\phi}{\phi}) \ln y   - \ln [\Gamma(1 / \phi)]
    \right ]


    &= \exp \left [
        \frac{y(1/\mu) - (-\ln \mu)}{-\phi}   - \frac{ \ln  \phi }{ \phi}  + (\frac{1-\phi}{\phi}) \ln y   - \ln [\Gamma(1 / \phi)]
    \right ]






和 :eq:`eq_gamma_011` 对比下，可以直接得到各个重要组件的形式。

.. math::

    \theta &= 1/\mu

    b(\theta) &= -\ln(\mu)

    a(\phi) &= - \phi

显然 Gamma模型的标准连接函数就是倒数函数， :math:`\eta=g(\mu) = \theta= 1/\mu` 。
现在我们看下 Gamma 分布的期望和方差函数。

.. math::

    b'(\theta) &= \frac{\partial b}{\partial \mu} \frac{\partial \mu}{\partial \theta}

    &= \left ( - \frac{1}{\mu} \right )(-\mu^2)

    &=\mu

    b''(\theta) &= \frac{\partial^2 b}{\partial \mu^2} \left ( \frac{\partial \mu}{\partial \theta} \right )
    + \frac{\partial b}{\partial \mu} \frac{\partial^2\mu}{\partial \theta^2}

    &= (1)(-\mu^2)

    &= -\mu^2


注意， :math:`b''(\theta)` 是方差函数，体现的是方差和均值的关系，
显然 **Gamma 分布的方差是和其均值相关的** ，Gamma 分布的方差为：

.. math::
    Var(y) = b''(\theta) a(\phi) = -\mu^2(-\phi) = \phi\mu^2

最后我们整理一下 Gamma 模型的一些关键组件。


.. math::

    \text{标准连接函数：}  & \eta = g(\mu) = \frac{1}{\mu}

    \text{反链接(响应)函数：}  & \mu = r(\eta) = \frac{1}{\eta}

    \text{方差函数：}  & \nu(\mu)=\mu^2

    \text{分散函数：}  & a(\phi) = -\phi

    \text{标准连接函数导数：} & g'=-\frac{1}{\mu^2}

参数估计
########################

似然函数
=========================

概率分布函数的指数族形式 :eq:`eq_gamma_013` ，直接去掉底数就得到了其对数似然函数。

.. math::

    \ell(\mu,\phi;y) = \sum_{i=1}^N \left \{
    \frac{y_i / \mu_i -(- \ln \mu_i)}{-\phi}
    + \frac{ 1-\phi}{\phi} \ln y_i
    - \frac{\ln \phi}{ \phi }
    - \ln \Gamma (1/ \phi)
    \right \}


根据 :eq:`eq_glm_estimate_ll_score`
，标准连接函数的Gamma模型的似然函数的一阶偏导为

.. math::

    U_j = \frac{\partial \ell}{\partial \beta_j}
    &= \sum_{i=1}^N \frac{y_i-\mu_i}{a(\phi) \nu(\mu_i)g(\mu_i)' }  x_{ij}

    &= - \sum_{i=1}^N \frac{y_i-\mu_i}{-\phi \mu_i^2 \eta_i^2 }  x_{ij}


IRLS
=========================

只需要给出 :math:`W` 和 :math:`Z` 的计算等式就可以应用IRLS算法。

.. math::

    W &= \text{diag} \left \{ \frac{ 1}{ a(\phi) \nu(\hat{\mu}) ( g' )^2}
    \right \}_{(N\times N)}

    &= \text{diag} \left \{ \frac{ -\hat{\mu}^2}{ \phi}
    \right \}_{(N\times N)}


.. math::


    Z &=   \left \{ (y- \hat{\mu}) g'  + \eta
    \right \}_{(N\times 1 )}

     &=   \left \{ \frac{-(y- \hat{\mu})}{ \hat{\mu}^2}  + \eta
    \right \}_{(N\times 1 )}


拟合优度
=========================

Gamma模型的偏差统计量为

.. math::

    D &= 2 \{ \ell(y;y)  - \ell(\hat{\mu};y)\}

    &= 2 \sum_{i=1}^N \frac{1+\ln(y_i)-y_i/\hat{\mu}_i-\ln(\hat{\mu}_i)  }{-\phi}

    &=   \frac{2}{\phi} \sum_{i=1}^N
     \left \{ \frac{y_i - \hat{\mu}_i}{\hat{\mu}_i} - \ln \left (  \frac{y_i}{\hat{\mu}_i} \right ) \right \}


Gamma模型的皮尔逊卡方统计量为


.. math::


    \chi^2
    &= \sum_{i=1}^N \frac{ (y_i-\hat{\mu}_i)^2}{\nu{\hat{\mu}_i}}

    &= \sum_{i=1}^N \frac{ (y_i-\hat{\mu}_i)^2}{\hat{\mu}_i^2}






其他连接函数
##########################################


对数 Gamma 模型
================================


前面我们提到，在给定一组特定的解释变量或预测变量的情况下，倒数链接估计模型响应的每单位速率。
对数链接的 Gamma 表示响应的对数率。
该模型规范与指数回归相同。
当然，这样的规范估计数据呈负指数下降。
但是，与生存分析中发现的指数模型不同，我们不能将对数伽马模型用于审查数据。
但是，我们看到未经审查的指数模型可以符合GLM规范。 我们将其保留到本章末尾。

对数伽玛模型，与它的对等模型一样，用于响应大于0的数据。几乎在每个学科中都可以找到示例。 例如，在健康分析中，通常可以使用对数伽玛回归来估算住院天数（LOS），
因为住院天数总是被约束为正数。
LOS数据通常使用泊松或负二项式回归进行估计，因为LOS的元素是离散的。
但是，当存在许多LOS元素（即许多不同的LOS值）时，许多研究人员发现伽马或高斯逆模型是可以接受的并且是更可取的。


在GLM之前，通常使用对数转换响应的高斯回归来估计现在使用对数伽马技术估计的数据。
尽管两种方法的结果通常相似，但对数伽马技术不需要外部转换，更易于解释，并带有一组残差，可用于评估模型的价值。
因此，对数伽马技术正在曾经使用过高斯技术的研究人员中得到越来越多的使用。


对数 gamma 模型是值其连接函数是对数函数，其响应函数(反链接)是指数函数。

.. math::

    \eta &= \ln (\mu)

    \mu &= e^{\eta}

对数连接函数的一阶导数是 :math:`g'(\mu)=1\mu` ，
可以轻松的使用IRLS算法进行参数估计。
但是由于对数连接函数表示规范连接函数，
IRLS算法的估计量和ML算法的估计量有不同的标准误差(standard errors)。
但是，除了极端情况外，标准误差的差异通常很小。
在较大数据集时，使用不同的估算方法通常不会产生任何推断差异。



恒等(identity) Gamma 模型
===================================================

恒等连接函数 :math:`\eta=\mu` ，假设 :math:`\mu` 和 :math:`\eta` 之间存在一一对应的关系。
高斯模型的规范链接就是恒等函数，但是Gamma模型中恒等函数不是规范链接。


在同一模型家族的不同链接之间进行选择有时可能很困难。
McCullagh和Nelder（1989）支持最小偏差的方法。
他们还检查残差以观察拟合的紧密度以及标准化残差本身的正态性和独立性。
在此示例中，我们有两个我们要比较的非规范链接：对数和恒等链接。
在其他因素都相同的情况下，最好选择偏差最小的模型。

我们也可以使用其他统计检验方法来评估链接之间的差异，
比如BIC和AIC，这些检验值较低的模型是更好一些的。
关于BIC，如果两个模型的BIC之间的绝对差小于2，
则两个模型的差异是比较小的。
2和6之间差值时，两个模型之间就有了一定的区别。
而6和10之间的差值，则说明两个模型有了明显的区别。
绝对差值大于10，那就肯定是BIC值较小的那个模型更好了。