Hurdle 与截断计数模型-未完纪

statsmodels 从 0.14 版本开始加入 Hurdle 与截断计数模型。

Hurdle 模型由两部分构成：一部分描述零值，另一部分描述大于零的计数分布。零值模型把“计数为零”与“计数大于零”视为二元结果；非零计数部分使用零截断的计数模型。

原文版本支持以泊松分布和负二项分布作为零值模型与计数模型；Logit、Probit、GLM-Binomial 等二元模型尚不能用作零值部分。泊松—泊松 Hurdle 模型的一个优点是：两部分参数相同时，标准泊松模型是它的特例，因此可以通过简单的 Wald 检验比较 Hurdle 模型与泊松模型。

这里实现的二元模型是一个在 1 处右删失的模型，也就是只能观察到 0 或 1。

假设各观测独立分布，观测之间没有相关性时，可以分别估计零值模型，以及针对零截断数据的计数模型。参数估计的协方差矩阵于是呈分块对角结构，对角块来自两个子模型。原文版本尚未实现联合估计。

删失与截断计数模型主要是为支持 Hurdle 模型而开发的。不过，左截断计数模型也有其他用途。这里的右删失模型本身没有单独的研究重点，因为它只支持二元观测，而二元观测也可以用 GLM-Binomial、Logit 或 Probit 建模。

Hurdle 模型使用单一的 HurdleCountModel 类，通过选项指定两个子模型的分布。截断模型对应的类是 TruncatedLFPoisson 和 TruncatedLFNegativeBinomialP；“LF”表示在固定、与各观测无关的截断点上进行左截断。

import numpy as np

from statsmodels.discrete.discrete_model import (
    Poisson,
)
from statsmodels.discrete.truncated_model import (
    HurdleCountModel,
)

模拟 Hurdle 模型

示例显式模拟一个泊松—泊松 Hurdle 模型，因为原文版本尚没有用于该模型的分布辅助函数。

np.random.seed(987456348)
# large sample to get strong results
nobs = 5000
x = np.column_stack((np.ones(nobs), np.linspace(0, 1, nobs)))

mu0 = np.exp(0.5 * 2 * x.sum(1))
y = np.random.poisson(mu0, size=nobs)
print(np.bincount(y))
y_ = y
indices = np.arange(len(y))
mask = mask0 = y > 0
for _ in range(10):

    print(mask.sum())
    indices = mask  # indices[mask]
    if not np.any(mask):
        break
    mu_ = np.exp(0.5 * x[indices].sum(1))
    y[indices] = y_ = np.random.poisson(mu_, size=len(mu_))
    np.place(y, mask, y_)
    mask = np.logical_and(mask0, y == 0)

np.bincount(y)

原文示例输出（本任务未执行）：

[102 335 590 770 816 739 573 402 265 176 116  59  35   7  11   4]
4898
602
93
11
2
0

原文示例输出（本任务未执行）：

array([ 102, 1448, 1502, 1049,  542,  234,   81,   31,    6,    5])

估计设定不正确的泊松模型

生成的数据具有零值偏少（zero deflation）的特征：观察到的零值数量少于泊松模型预期的数量。

拟合模型后，可以调用泊松诊断类的绘图函数，把预期的预测分布与实际频数进行比较。原图显示，泊松模型高估零值的数量，并低估计数为 1 和 2 的数量。

mod_p = Poisson(y, x)
res_p = mod_p.fit()
print(res_p.summary())

原文示例输出（本任务未执行）：

Optimization terminated successfully.
         Current function value: 1.668079
         Iterations 4
                          Poisson Regression Results                          
==============================================================================
Dep. Variable:                      y   No. Observations:                 5000
Model:                        Poisson   Df Residuals:                     4998
Method:                           MLE   Df Model:                            1
Date:                Thu, 27 Aug 2026   Pseudo R-squ.:                0.008678
Time:                        06:57:31   Log-Likelihood:                -8340.4
converged:                       True   LL-Null:                       -8413.4
Covariance Type:            nonrobust   LLR p-value:                 1.279e-33
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.6532      0.019     33.642      0.000       0.615       0.691
x1             0.3871      0.032     12.062      0.000       0.324       0.450
==============================================================================

dia_p = res_p.get_diagnostic()
dia_p.plot_probs();

泊松模型预测概率与实际计数频率对比

估计 Hurdle 模型

接下来，示例估计设定正确的泊松—泊松 Hurdle 模型。

HurdleCountModel 的函数签名与选项表明，泊松—泊松是默认组合，因此创建模型时无需另外指定分布选项：

HurdleCountModel(endog, exog, offset=None, dist='poisson', zerodist='poisson', p=2, pzero=2, exposure=None, missing='none', **kwargs)

HurdleCountModel 的结果类提供 get_diagnostic 方法，但原文版本只实现了部分诊断方法。预测分布图显示模型与数据高度一致。

mod_h = HurdleCountModel(y, x)
res_h = mod_h.fit(disp=False)
print(res_h.summary())

原文示例输出（本任务未执行）：

                     HurdleCountModel Regression Results                      
==============================================================================
Dep. Variable:                      y   No. Observations:                 5000
Model:               HurdleCountModel   Df Residuals:                     4996
Method:                           MLE   Df Model:                            2
Date:                Thu, 27 Aug 2026   Pseudo R-squ.:                 0.01503
Time:                        06:57:34   Log-Likelihood:                -8004.9
converged:               [True, True]   LL-Null:                       -8127.1
Covariance Type:            nonrobust   LLR p-value:                 8.901e-54
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
zm_const       0.9577      0.048     20.063      0.000       0.864       1.051
zm_x1          1.0576      0.121      8.737      0.000       0.820       1.295
const          0.5009      0.024     20.875      0.000       0.454       0.548
x1             0.4577      0.039     11.882      0.000       0.382       0.533
==============================================================================

dia_h = res_h.get_diagnostic()
dia_h.plot_probs();

Hurdle模型预测概率与实际计数频率对比

可以用 Wald 检验判断零值模型的参数是否与零截断计数模型的参数相同。原文示例中的 p 值很小，正确拒绝了“该模型只是泊松模型”的假设。这里使用大样本，所以该情形下检验的功效较高。

res_h.wald_test("zm_const = const, zm_x1 = x1", scalar=True)

原文示例输出（本任务未执行）：

<class 'statsmodels.stats.contrast.ContrastResults'>
<Wald test (chi2): statistic=470.6732075439199, p-value=6.231772522804731e-103, df_denom=2>

预测

Hurdle 模型可以预测整体模型及两个子模型的统计量，使用关键字 which 指定需要预测的统计量。

下面按原文 predict 文档字符串列出可用选项。which 是可选的字符串参数，默认值为 'mean'：

选项	含义
`'mean'`	因变量的条件期望 `E(y \| x)`。
`'mean-main'`	截断计数模型的均值参数。它不是截断后分布的均值。
`'linear'`	截断计数模型的线性预测值。
`'var'`	模型隐含的因变量估计方差。
`'prob-main'`	选择计数主模型的概率，也就是观察到非零计数的概率 `P(y > 0 \| x)`。
`'prob-zero'`	观察到零计数的概率 `P(y = 0 \| x)`，等于 `1 - prob-main`。
`'prob-trunc'`	截断计数模型的截断概率，也就是该截断模型隐含的零计数概率。
`'mean-nonzero'`	给定观测大于零时的条件期望 `E(y \| X, y > 0)`。
`'prob'`	从 0 到 `max(endog)` 的各计数概率；如果提供 `y_values`，则返回对应计数的概率。它返回多个值，对多个观测进行预测时结果为二维数组。

结果类的 predict 与 get_prediction 方法都支持这些选项。

接下来的示例，从原始数据中按等间隔选取解释变量，构成一组用于预测的解释变量，然后以这些变量为条件预测可用统计量。

which_options = [
    "mean",
    "mean-main",
    "linear",
    "mean-nonzero",
    "prob-zero",
    "prob-main",
    "prob-trunc",
    "var",
    "prob",
]
ex = x[slice(None, None, nobs // 5), :]
ex

原文示例输出（本任务未执行）：

array([[1.        , 0.        ],
       [1.        , 0.20004001],
       [1.        , 0.40008002],
       [1.        , 0.60012002],
       [1.        , 0.80016003]])

for w in which_options:
    print(w)
    pred = res_h.predict(ex, which=w)
    print("    ", pred)

原文示例输出（本任务未执行）：

mean
     [1.89150663 2.07648059 2.25555158 2.43319456 2.61673457]
mean-main
     [1.65015181 1.8083782  1.98177629 2.17180081 2.38004602]
linear
     [0.50086729 0.59243042 0.68399356 0.77555669 0.86711982]
mean-nonzero
     [2.04231955 2.16292424 2.29857565 2.45116551 2.62277411]
prob-zero
     [0.07384394 0.0399661  0.01871771 0.00733159 0.00230273]
prob-main
     [0.92615606 0.9600339  0.98128229 0.99266841 0.99769727]
prob-trunc
     [0.19202076 0.16391977 0.1378242  0.11397219 0.09254632]
var
     [1.43498239 1.51977118 1.63803729 1.7971727  1.99738345]
prob
     [[7.38439416e-02 3.63208532e-01 2.99674608e-01 1.64836199e-01
  6.80011882e-02 2.24424568e-02 6.17224344e-03 1.45501981e-03
  3.00125448e-04 5.50280612e-05]
 [3.99660987e-02 3.40376213e-01 3.07764462e-01 1.85518182e-01
  8.38717591e-02 3.03343722e-02 9.14266959e-03 2.36191491e-03
  5.33904431e-04 1.07277904e-04]
 [1.87177088e-02 3.10869602e-01 3.08037002e-01 2.03486809e-01
  1.00816333e-01 3.99590837e-02 1.31983274e-02 3.73659033e-03
  9.25635762e-04 2.03822556e-04]
 [7.33159258e-03 2.77316512e-01 3.01138113e-01 2.18003999e-01
  1.18365316e-01 5.14131777e-02 1.86098635e-02 5.77384524e-03
  1.56745522e-03 3.78244503e-04]
 [2.30272798e-03 2.42169151e-01 2.88186862e-01 2.28632665e-01
  1.36039066e-01 6.47558475e-02 2.56869828e-02 8.73374304e-03
  2.59833880e-03 6.87129546e-04]]

for w in which_options[:-1]:
    print(w)
    pred = res_h.get_prediction(ex, which=w)
    print("    ", pred.predicted)
    print("  se", pred.se)

原文示例输出（本任务未执行）：

mean
     [1.89150663 2.07648059 2.25555158 2.43319456 2.61673457]
  se [0.07877461 0.05693768 0.05866892 0.09551274 0.15359057]
mean-main
     [1.65015181 1.8083782  1.98177629 2.17180081 2.38004602]
  se [0.03959242 0.03164634 0.02471869 0.02415162 0.03453261]
linear
     [0.50086729 0.59243042 0.68399356 0.77555669 0.86711982]
  se [0.04773779 0.03148549 0.02960421 0.04397859 0.06453261]
mean-nonzero
     [2.04231955 2.16292424 2.29857565 2.45116551 2.62277411]
  se [0.02978486 0.02443098 0.01958745 0.0196433  0.02881753]
prob-zero
     [0.07384394 0.0399661  0.01871771 0.00733159 0.00230273]
  se [0.00918583 0.00405155 0.00220446 0.00158494 0.00090255]
prob-main
     [0.92615606 0.9600339  0.98128229 0.99266841 0.99769727]
  se [0.00918583 0.00405155 0.00220446 0.00158494 0.00090255]
prob-trunc
     [0.19202076 0.16391977 0.1378242  0.11397219 0.09254632]
  se [0.00760257 0.00518746 0.00340683 0.00275261 0.00319587]
var
     [1.43498239 1.51977118 1.63803729 1.7971727  1.99738345]
  se [0.04853902 0.03615054 0.02747485 0.02655145 0.03733328]

原文示例输出（本任务未执行）：

/opt/hostedtoolcache/Python/3.14.7/x64/lib/python3.14/site-packages/statsmodels/discrete/discrete_model.py:5286: UserWarning: using default log-link in get_prediction
  res = pred.get_prediction(

选项 which="prob" 为预测解释变量 exog 的每一行返回一个预测概率数组。实际分析中，经常关心对全部 exog 取平均后的概率。预测方法的 average=True 选项用于计算跨观测的预测值均值，以及这些平均预测的相应标准误与置信区间。

pred = res_h.get_prediction(ex, which="prob", average=True)
print("    ", pred.predicted)
print("  se", pred.se)

原文示例输出（本任务未执行）：

     [2.84324139e-02 3.06788002e-01 3.00960210e-01 2.00095571e-01
 1.01418732e-01 4.17809876e-02 1.45620174e-02 4.41222267e-03
 1.18509193e-03 2.86300514e-04]
  se [2.81472152e-03 5.00830805e-03 1.37524761e-03 1.87343644e-03
 1.99068656e-03 1.23878529e-03 5.78099178e-04 2.21180110e-04
 7.25021181e-05 2.08872555e-05]

示例用 pandas DataFrame 提供更易阅读的显示。predicted 列是响应值预测分布的概率质量函数，已对 exog 的 5 个网格点取平均。表中概率相加并不等于 1，因为比已观察计数更大的取值仍有正概率，而这些尾部取值没有列进表中；本例的尾部概率较小。

dfp_h = pred.summary_frame()
dfp_h

原文示例输出（本任务未执行）：

   predicted        se  ci_lower  ci_upper
0   0.028432  0.002815  0.022916  0.033949
1   0.306788  0.005008  0.296972  0.316604
2   0.300960  0.001375  0.298265  0.303656
3   0.200096  0.001873  0.196424  0.203767
4   0.101419  0.001991  0.097517  0.105320
5   0.041781  0.001239  0.039353  0.044209
6   0.014562  0.000578  0.013429  0.015695
7   0.004412  0.000221  0.003979  0.004846
8   0.001185  0.000073  0.001043  0.001327
9   0.000286  0.000021  0.000245  0.000327

prob_larger9 = pred.predicted.sum()
prob_larger9, 1 - prob_larger9

原文示例输出（本任务未执行）：

(np.float64(0.9999215487936677), np.float64(7.84512063323195e-05))

这里的 get_prediction 返回基类 PredictionResultsDelta 的一个实例。

对于依赖多个分布参数的非线性函数，标准误、p 值和置信区间等推断统计量通过 Delta 方法计算。预测推断以正态分布为基础。

pred

原文示例输出（本任务未执行）：

<statsmodels.base._prediction_inference.PredictionResultsDelta at 0x7fcdd9bd6470>

pred.dist, pred.dist_args

原文示例输出（本任务未执行）：

(<scipy.stats._continuous_distns.norm_gen at 0x7fcde0934440>, ())

可以把 Hurdle 模型预测的分布与之前估计的泊松模型预测分布进行比较。最后的 diff 列显示：在 exog 网格上的平均结果中，泊松模型对零值频数的高估约为观测数的 8%，对计数 1 和 2 的频数分别低估约 7% 与 3.7%。这些数字来自原文模拟示例。

pred_p = res_p.get_prediction(ex, which="prob", average=True)
dfp_p = pred_p.summary_frame()
dfp_h["poisson"] = dfp_p["predicted"]
dfp_h["diff"] = dfp_h["poisson"] - dfp_h["predicted"]
dfp_h

原文示例输出（本任务未执行）：

   predicted        se  ci_lower  ci_upper   poisson      diff
0   0.028432  0.002815  0.022916  0.033949  0.107848  0.079416
1   0.306788  0.005008  0.296972  0.316604  0.237020 -0.069768
2   0.300960  0.001375  0.298265  0.303656  0.263523 -0.037437
3   0.200096  0.001873  0.196424  0.203767  0.197657 -0.002439
4   0.101419  0.001991  0.097517  0.105320  0.112511  0.011093
5   0.041781  0.001239  0.039353  0.044209  0.051833  0.010052
6   0.014562  0.000578  0.013429  0.015695  0.020124  0.005561
7   0.004412  0.000221  0.003979  0.004846  0.006769  0.002356
8   0.001185  0.000073  0.001043  0.001327  0.002012  0.000827
9   0.000286  0.000021  0.000245  0.000327  0.000537  0.000250

其他估计后分析

估计后的 Hurdle 模型可以进行参数的 Wald 检验与预测，也提供对数似然值、信息准则等其他最大似然统计量。

原文版本尚未提供部分估计后方法。这些方法依赖的辅助函数并不是估计、参数推断与预测所必需的。主要尚不支持的方法为 score_test、get_distribution 和 get_influence。诊断度量仅支持以预测为基础的统计量。原文此处写作 get_diagnostics，前面的示例代码实际使用 get_diagnostic()；本文保留代码中的方法名。

res_h.llf, res_h.df_resid, res_h.aic, res_h.bic

原文示例输出（本任务未执行）：

(np.float64(-8004.904002793644),
 4996,
 np.float64(16017.808005587289),
 np.float64(16043.876778352953))

是否存在过度离散？可以使用 Pearson 残差计算按自由度归一化的 Pearson 卡方统计量；如果模型设定正确，这个比值应接近 1。

(res_h.resid_pearson**2).sum() / res_h.df_resid

原文示例输出（本任务未执行）：

np.float64(0.9989670114949286)

诊断类还提供诊断图使用的预测分布。原文版本目前没有提供其他统计量或检验。

dia_h.probs_predicted.mean(0)

原文示例输出（本任务未执行）：

array([0.02044612, 0.29147174, 0.29856288, 0.20740118, 0.10990976,
       0.04737579, 0.0172898 , 0.00548983, 0.00154646, 0.00039214])

res_h.resid[:10]

原文示例输出（本任务未执行）：

array([ 1.10849337,  1.10830496, -0.89188344, -0.89207183,  1.10773978,
       -0.8924486 , -0.89263697,  0.10717466,  0.1069863 ,  0.10679794])

正文、示例代码与原始图表来自官方仓库中的同名 notebook，按仓库 BSD 3-Clause 许可证使用。本文将全部正文单元译为中文，将预测选项改排为表格，并澄清 Pearson 统计量按自由度归一化；20 个代码单元与原文示例输出保留，图表从官方发布 notebook 提取。本任务未运行模型，原文结果的时间、路径和内存地址不表示本机实验。

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

a. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. b. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. c. Neither the name of statsmodels nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL STATSMODELS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

文章版权归作者所有，未经允许请勿转载。

THE END