使用 scipy.stats 将 Weibull 分布拟合到数据是不是表现不佳？

Posted 2023-02-19

【中文标题】使用 scipy.stats 将 Weibull 分布拟合到数据是不是表现不佳？

【英文标题】：Does fitting Weibull distribution to data using scipy.stats perform poor?使用 scipy.stats 将 Weibull 分布拟合到数据是否表现不佳？

【发布时间】：2021-01-08 19:33:46

【问题描述】：

我正在对一些整数数据拟合 Weibull 分布，并估计相关的形状、比例、位置参数。但是，我注意到 scipy.stats 库在这样做时性能不佳。

所以，我采取了不同的方向，并使用下面的代码检查了合身性能。我首先使用 Weibull 分布创建 100 个数字，参数 shape=3、scale=200、location=1。随后，我使用 fitter 库估计最佳分布拟合。

from fitter import Fitter
import numpy as np
from scipy.stats import weibull_min

# generate numbers
x = weibull_min.rvs(3, scale=200, loc=1, size=100)

# make them integers
data = np.asarray(x, dtype=int)

# fit one of the four distributions
f = Fitter(data, distributions=["gamma", "rayleigh", "uniform", "weibull_min"])
f.fit()

f.summary()

我希望最适合的是 Weibull 分布。我试过重新运行这个测试。有时 Weibull 拟合是一个很好的估计。然而，大多数时候 Weibull 拟合被报告为最差的结果。在这种情况下，估计的参数为 = (0.13836651040093312, 66.99999999999999, 1.3200752378443505)。我假设这些参数依次对应于形状、比例、位置。以下是拟合过程的摘要。

$ f.summary()
             sumsquare_error          aic          bic  kl_div
gamma               0.001601  1182.739756 -1090.410631     inf
rayleigh            0.001819  1154.204133 -1082.276256     inf
uniform             0.002241  1113.815217 -1061.400668     inf
weibull_min         0.004992  1558.203041  -976.698452     inf

另外，生成了以下情节。

此外，瑞利分布是形状参数 = 2 的 Weibull 分布的特例。因此，我希望得到的 Weibull 拟合至少与 Rayleigh 一样好。

更新

我在 numpy 版本 1.19.2 和 scipy 版本 1.5.2 的 Linux/Ubuntu 20.04 机器上运行了上述测试。上面的代码似乎按预期运行，并为 Mac 机器上的 Weibull 分发返回了正确的结果。

我还测试了使用 R 库 fitdistrplus 在 Linux 机器上生成的数据 x 拟合 Weibull 分布：

fit.weib <- fitdist(x, "weibull")

并观察到估计的形状和比例值与最初给定的值非常接近。到目前为止，最好的猜测是问题是由于一些 Python-Ubuntu 错误/不兼容造成的。

我可以被认为是这个领域的新手。所以，我想知道，我在这里做错了吗？还是以某种方式预期这个结果？非常感谢任何帮助。

谢谢。

【参考方案1】：

库fitter 不允许为诸如 a、loc 等发行版指定参数。奇怪的是，Mac 产生了更好的拟合，而 Linux 则为相同版本的 Numpy 和 Scipy 提供了最合适的结果。根本原因可能包括为 Linux 和 Mac 设计的不同 BLAS-LAPACK 算法，https://***.com/a/49274049/6806531 或 weibull_min 可能无法初始化在线讨论的参数a = 1，或默认浮点精度。但是，可以解决fitter 库中的错误。知道weib_min是expon_weib，参数a固定为1，把fitter.py中_timed_run函数里面的run函数改成

        def run(self):

            try:
                if distribution == "exponweib":
                    self.result = func(args,floc=0,fa = 1, **kwargs)
                else:
                    self.result = func(args, floc=0, **kwargs)
            except Exception as err:
                self.exc_info = sys.exc_info()

使用 exponweib 作为 weib_min 的结果与 R fitdist 几乎相同。

【讨论】：

谢谢。在您指出的 if-else 块中，使用以下任一方法： if distribution == "exponweib": self.result = func(args,floc=0,fa = 1, **kwargs) 或 if distribution == "weibull_min ": self.result = func(args, floc=0, **kwargs) 我得到了预期的估计参数值。

【参考方案2】：

我不熟悉 Fitter 库，但为了得出一些结论，我建议：

重试您的代码，但采用 size=10,000。在这种情况下，有足够的数据点供拟合方法使用。从理论上讲，您会期望 Weibull 提供最合适的结果。

我注意到位置参数有时会很痛苦。您可以尝试通过使用 floc=1 固定位置参数来运行您的拟合（即等于您的位置采样参数）。你得到了什么？此外，仅供参考，使用 MLE，只需 loc=min(x)，其中 x 是您的数据集。对于指数分布，这实际上是位置参数的 MLE。对于其他发行版我不确定，但如果这也适用于其他发行版，我不会感到惊讶。这将减少 1 个参数的拟合过程。

最后，我注意到，如果您对某些分布的位置/比例/形状取较小的值，则 scipy.stats 分布的函数 logpdf 和 logcdf 会产生 np.inf 值。在这种情况下，您或许可以使用 Powell 优化算法并设置参数值的界限。