计算性能 scipy weibull min fit vs Matlab wblfit

Posted 2023-02-19

【中文标题】计算性能 scipy weibull min fit vs Matlab wblfit

【英文标题】：Calculation performance scipy weibull min fit vs Matlab wblfit

【发布时间】：2021-10-16 11:45:47

【问题描述】：

使用 Matlab wblrnd 和 wblfit 函数以及 Python scipy.stats.weibull_min.fit 函数将数据拟合到 Weibull 分布，我发现 Matlab 的性能比 Python 高出近 2 个数量级。我正在寻求一些帮助来提高 Python 代码的性能。

问题：

在将 Matlab 代码转换为 Python 时，我遇到了以下代码：

weibull_parameters = zeros(10000, 2)
for i = 1:10000
    data = sort(wblrnd(alpha, beta, 1, 24))
    [weibull_parameters(i, :), ~] = wblfit(data, confidence_interval, censoring_array)
end

此代码从 Weibull 分布生成 24 个随机数，然后将生成的数据向量再次拟合到 Weibull 分布。

在 Python 中，我将其翻译为：

from scipy.stats import weibull_min
import numpy as np


data = np.sort(alpha * np.random.default_rng().weibull(beta, (10000, 24)))

weibull_parameters = np.zeros((10000, 2))
for idx, row in enumerate(data):
    weibull_parameters[idx, :] = weibull_min.fit(row, floc=0)[2::-2]

在这里，我一次生成完整的随机数据，然后使用 weibull_min.fit 函数遍历行以获取相应的 Weibull 参数。最后的切片就是在输出中只选择尺度和形状参数，并按正确的顺序排列。

我遇到的主要问题是Python中的计算性能很糟糕。 Matlab 在几秒钟内运行此代码，但对于 Python，每 100 次迭代需要 1-1.5 秒（在我的笔记本电脑上），因此性能差异几乎是 2 个数量级。

有什么方法可以提高 Python 的性能吗？是否可以矢量化拟合计算？很遗憾，我在网上找不到关于此主题的任何内容。

注意 1：Matlab 允许用户在 wblfit 函数中指定置信区间，但是对于 Python，我找不到包含它的方法，所以我忽略了它。

注意 2：我能找到的唯一包含审查的选项是使用 surpyval package，但性能更加糟糕（每 100 次迭代大约需要 10 秒）

【参考方案1】：

Python 并不是最快的语言。您可以采取一些措施来加快速度，但您会发现准确性和速度之间存在平衡。

至于适合 Weibull 分布的方法，有几个包可以做到这一点。包scipy、surpyval、lifelines 和reliability 都适合完整的数据。最后 3 个也将处理 scipy 不会处理的审查数据。

我是可靠性的作者，所以我将向您展示一个使用此包的示例：

from reliability.Distributions import Weibull_Distribution
from reliability.Fitters import Fit_Weibull_2P
import time
import numpy as np

rows=100
samples = 24
data_array = np.empty((rows,samples))
true_parameters = np.empty((rows,2))
for i in range(rows):
    alpha = np.random.randint(low=1,high=999) + np.random.rand() #alpha between 1 and 1000
    beta = np.random.randint(low=1,high=10) - np.random.rand()/2 #beta between 0.5 and 10
    true_parameters[i][0] = alpha
    true_parameters[i][1] = beta
    dist = Weibull_Distribution(alpha=alpha,beta=beta)
    data_array[i]=dist.random_samples(samples)

start_time = time.time()
parameters = np.empty((rows,2))
for i in range(rows):
    fit = Fit_Weibull_2P(failures=data_array[i],show_probability_plot=False,print_results=False)
    parameters[i][0] = fit.alpha
    parameters[i][1] = fit.beta

runtime = time.time()-start_time
# np.set_printoptions(suppress=True) #supresses the scientific notation used by numpy
# print('True parameters:')
# print(true_parameters)
# print('Fitted parameters:')
# print(parameters)
print('Runtime:',runtime,'seconds')
print('Runtime per iteration:',runtime/rows,'seconds')

当我运行它时，它会给出：

Runtime: 3.378781318664551 seconds
Runtime per iteration: 0.033787813186645504 seconds

根据您在问题中引用的时间，这大约是 scipy 的两倍，但只有 surpyval 所用时间的三分之一。

我希望这有助于向您展示一种不同的方式来做同样的事情，但我知道这可能不是您所寻求的性能改进。获得大幅性能提升的唯一方法是在纯 python 中使用最小二乘估计，也许使用 numba 加速。这种方法可能会给您带来不如 MLE 的结果，但正如我之前所说，速度和准确性之间以及速度和编码便利性之间存在平衡。

【讨论】：

感谢您向我展示可靠性包，我确实不知道。不幸的是，我不能放弃使用 MLE，因为这是行业标准。我仍然对 MATLAB 能够如此快速地执行拟合感到震惊。也许它在后台做了很多并行化/矢量化。现在，虽然您显示的结果表明至少有一种更有利的实施审查方法，但性能仍然太慢。因为我们需要在代码中运行拟合部分 4 次，所以这意味着人们要等待 20 分钟才能完成代码。也许 Numba 有帮助

如果您需要使用 MLE，那么 numba 将无法帮助您（除非您愿意深入研究 autograd 的源代码）。您的数据集是否包含大量重复值？如果是这样，那么您可能会更快地找到生命线。