如何使用 Python 获得 Weibull 分布的置信区间？

Posted 2023-02-19

【中文标题】如何使用 Python 获得 Weibull 分布的置信区间？

【英文标题】：How to get the confidence interval of a Weibull distribution using Python?

【发布时间】：2019-12-26 16:12:59

【问题描述】：

我想通过 Python 执行具有 0.95% 置信区间的概率 Weibull 拟合。作为测试数据，我使用了针对可靠性 R(t) 绘制的测量失败周期。

到目前为止，我找到了一种执行 Weibull 拟合的方法，但是，我仍然无法获得置信界限。具有相同测试数据集的 Weibull 图已经使用原点执行，因此我知道我会“预期”置信区间的形状。但我不明白如何到达那里。

我在reliawiki 上找到了有关 Weibull 置信区间的信息（参见基于 Fisher 矩阵置信界限的可靠性界限），并使用那里的描述来计算方差和置信上限和下限（R_U 和 R_L）。

这是我的 Weibull 拟合的工作代码示例，以及我对基于 reliawiki 描述的测试数据集的置信界限（参见可靠性界限）。对于拟合，我使用了 OLS 模型拟合。

import os, sys
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy.optimize import curve_fit
import math
import statsmodels.api as sm

def weibull_ticks(y, pos):
    return ":.0f%".format(100 * (1 - np.exp(-np.exp(y))))

def loglog(x):
    return np.log(-np.log(1 - np.asarray(x)))

class weibull_example(object):

    def __init__(self, dat):
        self.fits = 
        dat.index = np.arange(1, len(dat) + 1)
        dat.sort_values('data', inplace=True)
        #define yaxis-values
        dat['percentile'] = dat.index*1/len(dat)
        self.data = dat

        self.fit()
        self.plot_data()

    def fit(self):
        #fit the data points with a the OLS model
        self.data=self.data[:-1]
        x0 = np.log(self.data.dropna()['data'].values)
        Y = loglog(self.data.dropna()['percentile'])
        Yx = sm.add_constant(Y)
        model = sm.OLS(x0, Yx)
        results = model.fit()
        yy = loglog(np.linspace(.001, .999, 100))
        YY = sm.add_constant(yy)
        XX = np.exp(results.predict(YY))
        self.eta = np.exp(results.params[0])
        self.beta = 1 / results.params[1]
        self.fits['syx'] = 'results': results, 'model': model,
                            'line': np.row_stack([XX, yy]),
                            'beta': self.beta,
                            'eta': self.eta

        cov = results.cov_params()
        #get variance and covariance
        self.beta_var = cov[1, 1]
        self.eta_var = cov[0, 0]
        self.cov = cov[1, 0]

    def plot_data(self, fit='yx'):
        dat = self.data
        #plot data points
        plt.semilogx(dat['data'], loglog(dat['percentile']), 'o')
        fit = 's' + fit
        self.plot_fit(fit)

        ax = plt.gca()
        formatter = mpl.ticker.FuncFormatter(weibull_ticks)
        ax.yaxis.set_major_formatter(formatter)
        yt_F = np.array([0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,
                         0.6, 0.7, 0.8, 0.9, 0.95, 0.99])
        yt_lnF = loglog(yt_F)
        plt.yticks(yt_lnF)

        plt.ylim(loglog([.01, .99]))

    def plot_fit(self, fit='syx'):
        dat = self.fits[fit]['line']
        plt.plot(dat[0], dat[1])

        #calculate variance to get confidence bound
        def variance(x):
            return (math.log(x) - math.log(self.eta)) ** 2 * self.beta_var + \
                   (self.beta/self.eta) ** 2 * self.eta_var - \
                   2 * (math.log(x) - math.log(self.eta)) * (-self.beta/self.eta) * self.cov

        #calculate confidence bounds
        def confidence_upper(x):
            return 1-np.exp(-np.exp(self.beta*(math.log(x)-math.log(self.eta)) - 0.95*np.sqrt(variance(x))))
        def confidence_lower(x):
            return 1-np.exp(-np.exp(self.beta*(math.log(x)-math.log(self.eta)) + 0.95*np.sqrt(variance(x))))

        yvals_1 = list(map(confidence_upper, dat[0]))
        yvals_2 = list(map(confidence_lower, dat[0]))

        #plot confidence bounds
        plt.semilogx(dat[0], loglog(yvals_1), linestyle="solid", color="black", linewidth=2,
                 label="fit_u_1", alpha=0.8)
        plt.semilogx(dat[0], loglog(yvals_2), linestyle="solid", color="green", linewidth=2,
                 label="fit_u_1", alpha=0.8)

def main():
    fig, ax1 = plt.subplots()
    ax1.set_xlabel("$Cycles\ til\ Failure$")
    ax1.set_ylabel("$Weibull\ Percentile$")

    #my data points
    data = pd.DataFrame('data': [1556, 2595, 11531, 38079, 46046, 57357])
    weibull_example(data)

    plt.savefig("Weibull.png")
    plt.close(fig)

if __name__ == "__main__":
    main()

我的情节中的置信区间看起来不像我预期的那样。我尝试了很多不同的“差异”，只是为了了解功能并检查问题是否只是打字错误。同时，我确信这个问题更普遍，我从reliawiki 的描述中理解了一些错误的东西。不幸的是，我真的不明白问题出在哪里，我不知道我可以问其他任何人。在互联网和不同的论坛上，我没有找到合适的答案。

这就是我决定在这里问这个问题的原因。这是我第一次在论坛上提问。因此，我希望我对所有内容都进行了充分的解释，并且代码示例很有用。 非常感谢:)

【参考方案1】：

对于迟到的答案深表歉意，但我会为任何未来的读者提供它。 与其尝试自己实现这一点，不如考虑使用专门为此而设计的包，称为可靠性。 这是您的用例的example。

如果对你有帮助，记得点赞这个答案:)