scipy.stats.weibull_min.fit() - 如何处理右删失数据？

Posted 2023-02-19

【中文标题】scipy.stats.weibull_min.fit() - 如何处理右删失数据？

【英文标题】：scipy.stats.weibull_min.fit() - how to deal with right-censored data?

【发布时间】：2018-05-28 02:34:56

【问题描述】：

非删失（完整）数据集

我正在尝试使用scipy.stats.weibull_min.fit() 函数来拟合一些生命数据。生成的示例数据包含在下面的 values 中。

values = np.array(
    [10197.8, 3349.0, 15318.6, 142.6, 20683.2, 
    6976.5, 2590.7, 11351.7, 10177.0, 3738.4]
)

我尝试使用该函数进行拟合：

fit = scipy.stats.weibull_min.fit(values, loc=0)

结果：

(1.3392877335100251, -277.75467055900197, 9443.6312323849124)

这与标称 beta 和 eta 值 1.4 和 10000 相差不远。

右删失数据

weibull 分布以其处理右删失数据的能力而闻名。这使得它对可靠性分析非常有用。如何处理scipy.stats 中的右删失数据？即，对尚未出现故障的数据进行曲线拟合？

输入表单可能如下所示：

values = np.array(
    [10197.8, 3349.0, 15318.6, 142.6, np.inf, 
    6976.5, 2590.7, 11351.7, 10177.0, 3738.4]
)

或者可能使用np.nan 或简单地使用0。

np 的两个解决方案都在抛出 RunTimeWarnings 并且绝对没有接近正确的值。我使用数值——比如0和-1——去掉了RunTimeWarning，但是返回的参数明显有缺陷。

其他软件

在一些可靠性或寿命分析软件（minitab，lifelines）中，需要有两列数据，一列是实际数字，一列表示该项目是否失败。例如：

values = np.array(
    [10197.8, 3349.0, 15318.6, 142.6, 0, 
    6976.5, 2590.7, 11351.7, 10177.0, 3738.4]
)

censored = np.array(
    [True, True, True, True, False,
    True, True, True, True, True]
)

我在文档中看不到这样的路径。

【问题讨论】：

在路上的某个地方，似乎使用了np.log，我猜inf、nan 和0 会导致问题出现。您是否可以以某种方式替换这些条目，values[np.isinf(values)] = 10000.、values[np.isnan(values)] = 1. 和 values[np.isclose(values, 0.)] = 10 ** (-6) 或类似的东西？

@Cleb 这不是数值问题。例如 0 在算法通过时工作得很好，但 0 根本不会产生正确的结果。

好吧，只是一个疯狂的猜测，正如我看到的RuntimeWarning: invalid value encountered in subtract return np.log(c) + sc.xlogy(c - 1, x) - pow(x, c)；所以我认为log(c) 部分可能是问题所在。然后，当我运行您的示例并将 np.inf 替换为另一个更高的值时，它运行良好。但我对此了解的太少，恐怕无法提供帮助……

正如您提到的“其他软件”：您能否更具体地说明一下，如果可以的话，还请提供您在其他语言中如何称呼它的代码？

@Cleb 其他软件有minitab之类的（我不是用户，但看过视频教程）。 lifelines 项目支持类似的接口。一段时间以来，我一直在为这种类型的分析尝试不同的途径，所以我也有一个github repository，我正在使用它来获得我需要的功能。我可能最终会尝试将其与生命线合并，但我希望在此之前让它更实用。

【参考方案1】：

老问题，但如果有人遇到这个问题，有一个新的 python 生存分析包 surpyval 可以处理这个问题以及其他审查和截断情况。对于您在上面提供的示例，它只是：

import surpyval as surv
values = np.array([10197.8, 3349.0, 15318.6, 142.6, 6976.5, 2590.7, 11351.7, 10177.0, 3738.4])

# 0 = failed, 1 = right censored
censored = np.array([0, 0, 0, 0, 0, 1, 1, 1, 0])

model = surv.Weibull.fit(values, c=censored)
print(model.params)

(10584.005910580288, 1.038163987652635)

您可能还对 Weibull 图感兴趣：

model.plot(plot_bounds=False)

Weibull plot

完全披露，我是surpyval的创造者