GenericLikelihoodModel 拟合二项分布：输出 stderr 为 nan答案

【问题标题】：GenericLikelihoodModel to fit binomial distribution: output stderr is nanGenericLikelihoodModel 拟合二项分布：输出 stderr 为 nan
【发布时间】：2017-12-31 18:46:50
【问题描述】：

由于 python 中没有直接的离散分布拟合包，我尝试使用 statsmodels.base.model.GenericLikelihoodModel 来拟合二项分布。但是，在某些情况下，例如 n=5，p=0.5，估计摘要中的每个项都是 nan，除了点估计。

如果我将 n 从 5 更改为 10，那么一切正常。

我的代码是：

import numpy as np
from statsmodels.base.model import GenericLikelihoodModel
from scipy.stats import binom


class Binom(GenericLikelihoodModel):
    def loglike(self, params):            
        return np.log(binom.pmf(self.endog, *params)).sum()

n, p = 5, 0.5
params = (n, p)
x = binom.rvs(5, p, size=1000, random_state=1)

res = Binom(x).fit(start_params=params)
res.df_model = len(params)
res.df_resid = len(x) - len(params)
print(res.summary())

相关输出：

        coef    std err          z      P>|z|      [95.0% Conf. Int.]
par0  5.0000        nan        nan        nan       nan       nan
par1  0.4972        nan        nan        nan       nan       nan

错误信息：

//anaconda/lib/python3.6/site-packages/ipykernel/main.py:23: RuntimeWarning: 除以零在日志中遇到

//anaconda/lib/python3.6/site-packages/statsmodels/tools/numdiff.py:329：RuntimeWarning：double_scalars中遇到无效值 - f(*((x - ee[i,:] - ee[j,:],)+args), **kwargs),)

//anaconda/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:875：RuntimeWarning：在更大时遇到无效值返回 (self.a

//anaconda/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:875：RuntimeWarning：在less中遇到无效值返回 (self.a

//anaconda/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:1814：RuntimeWarning：在less_equal中遇到无效值 cond2 = cond0 & (x

我终于找到了我的代码失败的原因。在statsmodels中，std err的计算是通过公式

np.sqrt(np.diag(self.cov_params()))

其中 cov_params 是点估计处的 hessian 矩阵。在拟合过程中，hessian矩阵的计算公式为

hess[i, j] = (f(*((x + ee[i, :] + ee[j, :],) + args), **kwargs)
                      - f(*((x + ee[i, :] - ee[j, :],) + args), **kwargs)
                      - (f(*((x - ee[i, :] + ee[j, :],) + args), **kwargs)
                      - f(*((x - ee[i, :] - ee[j, :],) + args), **kwargs),)
                      )/(4.*hess[i, j])

其中 f 是 loglike 函数，x 是点估计，ee 是估计误差。如果n的点估计是真实的n，那么需要为hess[i,j]计算f(n-ee, p)。 n-ee

结合我自定义的loglike函数，我们可以看到它涉及log(0)=Inf的计算（解释如下），而hess[i,j]涉及log(0)的加法（导致nan）。这就是为什么最终的 std err 和所有其他指标都变为 nan。

log(0) 中的零是通过 binom.pmf(x, n, p)=0 得到的，其中 x>n。在我的实践中，样本大小为 1000。鉴于 binom.pmf(5,5,0.5)=0.03125 和 binom.pmf(10,10,0.5)=0.00097，我们可以看到当我选择 n=10,p=0.5对于抽样，样本中可能不存在 10。因此我们不会遇到 log(0) 的情况。但是，如果我将样本大小增加到 1000000，则会出现预期的样本错误。

所以我现在的问题变成了有没有办法通过修改我的自定义 loglike 函数或通过覆盖 statsmodels 中的 hessian 函数来计算这种情况下的 std err？

【问题讨论】：

标签： python statsmodels mle

【解决方案1】：

一般来说，当分布的支持取决于参数时，最大似然估计存在理论和计算问题。

在这种情况下，n 定义了分布的上限，据我所知，标准误差等推理统计数据在这种情况下无效。

在通常的二项式情况下，我们假设试验次数 n 是已知的，我们只估计概率或比例。在 statsmodels 和其他软件包中，这可以通过使用带有家族二项式的 GLM 并将成功和失败的数量指定为因变量来估计。默认情况下，使用 Logit 链接，该链接要求对参数进行转换以获得概率，在这种情况下通过 predict 方法完成。

【讨论】：