1

假设我有一个样本，我有理由相信它遵循指数分布。我想估计分布参数 (lambda) 和一些置信度指示。置信区间或标准误差都可以。可悲的是，scipy.stats.expon.fit似乎不允许这样做。这是一个例子，我将使用 lambda=1/120=0.008333 作为测试数据：

""" Generate test data"""
import scipy.stats
test_data = scipy.stats.expon(scale=120).rvs(size=3000)

""" Scipy.stats fit"""
fit = scipy.stats.expon.fit(test_data)
print("\nExponential parameters:", fit, " Specifically lambda: ", 1/fit[1], "\n")

# Exponential parameters: (0.0066790678905608875, 116.8376079908356)  Specifically lambda:  0.008558887991599736

关于伽马分布数据的类似问题的答案建议使用GenericLikelihoodModel模块statsmodels。虽然我可以确认这对于 gamma 分布的数据非常有效，但它不适用于指数分布，因为优化显然会导致不可逆的 Hessian 矩阵。这是由于 Hessian 矩阵中的非有限元素或np.linalg.eigh为 Hessian 矩阵生成非正特征值造成的。（这里的源代码；HessianInversionWarning 在类的fit方法中引发LikelihoodModel。）。

""" Statsmodel fit"""
from statsmodels.base.model import GenericLikelihoodModel

class Expon(GenericLikelihoodModel):
    nparams = 2
    def loglike(self, params):
        return scipy.stats.expon.logpdf(self.endog, *params).sum()


res = Expon(test_data).fit(start_params=fit)
res.df_model = len(fit)
res.df_resid = len(test_data) - len(fit)
print(res.summary())

#Optimization terminated successfully.
#         Current function value: 5.760785
#         Iterations: 38
#         Function evaluations: 76
#/usr/lib/python3.8/site-packages/statsmodels/tools/numdiff.py:352: RuntimeWarning: invalid value encountered in double_scalars
#  hess[i, j] = (f(*((x + ee[i, :] + ee[j, :],) + args), **kwargs)
#/usr/lib/python3.8/site-packages/statsmodels/base/model.py:547: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available
#  warn('Inverting hessian failed, no bse or cov_params '
#/usr/lib/python3.8/site-packages/scipy/stats/_distn_infrastructure.py:903: RuntimeWarning: invalid value encountered in greater
#  return (a < x) & (x < b)
#/usr/lib/python3.8/site-packages/scipy/stats/_distn_infrastructure.py:903: RuntimeWarning: invalid value encountered in less
#  return (a < x) & (x < b)
#/usr/lib/python3.8/site-packages/scipy/stats/_distn_infrastructure.py:1912: RuntimeWarning: invalid value encountered in less_equal
#  cond2 = cond0 & (x <= _a)
#                                Expon Results                                 
#==============================================================================
#Dep. Variable:                      y   Log-Likelihood:                -17282.
#Model:                          Expon   AIC:                         3.457e+04
#Method:            Maximum Likelihood   BIC:                         3.459e+04
#Date:                Thu, 06 Aug 2020                                         
#Time:                        13:55:24                                         
#No. Observations:                3000                                         
#Df Residuals:                    2998                                         
#Df Model:                           2                                         
#==============================================================================
#                 coef    std err          z      P>|z|      [0.025      0.975]
#------------------------------------------------------------------------------
#par0           0.0067        nan        nan        nan         nan         nan
#par1         116.8376        nan        nan        nan         nan         nan
#==============================================================================

这似乎每次都会发生，因此它可能与指数分布的数据有关。

还有其他可能的方法吗？或者我可能在这里遗漏了什么或做错了什么？

编辑：事实证明，我做错了什么，即我错误地

test_data = scipy.stats.expon(120).rvs(size=3000)

代替

test_data = scipy.stats.expon(scale=120).rvs(size=3000)

并相应地查看 fit 元组的第一个元素，而我应该查看第二个元素。

结果，我考虑的另外两个选项（按照维基百科上描述的标准程序手动计算拟合和置信区间）并按照此答案scikits.bootstrap中的建议使用实际上确实有效，并且是我将在一分钟内添加的解决方案的一部分不是问题。

4

2 回答 2

1

如已编辑问题中所述，部分问题是我在创建样本时查看了错误的参数，并再次进行拟合。

剩下的是，scipy.stats.expon.fit它不提供计算置信度或错误的可能性，并且由于 Hessian 格式错误，使用此处建议GenericLikelihoodModel的statsmodels模块失败。

但是，有三种方法确实有效：

1. 使用维基百科文章中给出的指数数据置信区间的简单推理程序

""" Maximum likelihood"""
import numpy as np
ML_lambda = 1 / np.mean(test_data)
print("\nML lambda: {0:8f}".format(ML_lambda))

#ML lambda: 0.008558

""" Bias corrected ML"""
ML_BC_lambda = ML_lambda - ML_lambda / (len(test_data) - 1)
print("\nML bias-corrected lambda: {0:8f}".format(ML_BC_lambda))

#ML bias-corrected lambda: 0.008556

置信区间的计算：

""" Maximum likelihood 95% confidence"""
CI_distance = ML_BC_lambda * 1.96/(len(test_data)**0.5)
print("\nLambda with confidence intervals: {0:8f} +/- {1:8f}".format(ML_BC_lambda, CI_distance))
print("Confidence intervals: ({0:8f}, {1:9f})".format(ML_BC_lambda - CI_distance, ML_BC_lambda + CI_distance))

#Lambda with confidence intervals: 0.008556 +/- 0.000306
#Confidence intervals: (0.008249,  0.008862)

第二种选择：此外，置信区间方程也应该适用于由不同产生的 lambda 估计，例如来自 的scipy.stats.expon.fit。（我认为拟合过程scipy.stats.expon.fit更可靠，但事实证明它实际上是相同的，没有偏差校正（见上文）。）

""" Maximum likelihood 95% confidence based on scipy.stats fit"""
scipy_stats_lambda = 1 / fit[1]
scipy_stats_CI_distance = scipy_stats_lambda * 1.96/(len(test_data)**0.5)
print("\nOr, based on scipy.stats fit:")
print("Lambda with confidence intervals: {0:8f} +/- {1:8f}".format(scipy_stats_lambda, scipy_stats_CI_distance))
print("Confidence intervals: ({0:8f}, {1:9f})".format(scipy_stats_lambda - scipy_stats_CI_distance, 
                                                                scipy_stats_lambda + scipy_stats_CI_distance))

#Or, based on scipy.stats fit:
#Lambda with confidence intervals: 0.008559 +/- 0.000306
#Confidence intervals: (0.008253,  0.008865)

scikits.bootstrap2.遵循此答案中的建议进行引导 这会产生一个InstabilityWarning: Some values were NaN; results are probably unstable (all values were probably equal)，因此应该对此持怀疑态度。

""" Bootstrapping with scikits"""
print("\n")
import scikits.bootstrap as boot
bootstrap_result = boot.ci(test_data, scipy.stats.expon.fit)
print(bootstrap_result)

#tmp/expon_fit_test.py:53: InstabilityWarning: Some values were NaN; results are probably unstable (all values were probably equal)
#  bootstrap_result = boot.ci(test_data, scipy.stats.expon.fit)
#[[6.67906789e-03 1.12615588e+02]
# [6.67906789e-03 1.21127091e+02]]

3.使用rpy2

""" Using r modules with rpy2"""
import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
MASS = importr('MASS')
import rpy2.robjects.numpy2ri
rpy2.robjects.numpy2ri.activate()
rpy_fit = MASS.fitdistr(test_data, "exponential")
rpy_estimate = rpy_fit.rx("estimate")[0][0]
rpy_sd = rpy_fit.rx("sd")[0][0]
rpy_lower = rpy_estimate - 2*rpy_sd
rpy_upper = rpy_estimate + 2*rpy_sd
print("\nrpy2 fit: \nLambda={0:8f} +/- {1:8f}, CI: ({2:8f}, {3:8f})".format(rpy_estimate, rpy_sd, rpy_lower, rpy_upper))

#rpy2 fit: 
#Lambda=0.008558 +/- 0.000156, CI: (0.008246, 0.008871)
于 2020-08-06T13:41:39.613 回答
0

您已经找到了解决问题的方法，但这里有一个基于 OpenTURNS 的解决方案。我认为它依赖于引擎盖下的引导。

OpenTURNS 要求您重新塑造数据，以便清楚地看到我们正在处理 3000 个一维点，而不是单个 3000 维点。

test_data = test_data.reshape(-1, 1)

其余的比较简单。

import openturns as ot

confidence_level = 0.9

params_ci = ot.ExponentialFactory().buildEstimator(test_data).getParameterDistribution().computeBilateralConfidenceInterval(confidence_level)

lambda_ci = [params_ci.getLowerBound()[0], params_ci.getUpperBound()[0]] 
# the index 0 means we are interested in the CI on lambda

print(lambda_ci)

我得到以下输出（但这取决于随机种子）：

[0.008076302149561718, 0.008688296487447742]
于 2020-10-03T20:10:09.910 回答