Как получить доверительный интервал распределения Вейбулла с помощью Python?

Я хочу выполнить оценку вероятности Вейбулла с доверительным интервалом 0,95% с помощью Python. В качестве тестовых данных я использую неудачные циклы измерения, которые нанесены на график в зависимости от надежности R(t).

До сих пор я нашел способ выполнить подгонку Вейбулла, однако мне все еще не удается получить доверительные границы. График Вейбулла с тем же набором тестовых данных уже был выполнен с исходной точкой, поэтому я знаю, какую форму я «ожидаю» для доверительного интервала. Но я не понимаю, как туда попасть.

Я нашел информацию о доверительных интервалах Вейбулла на reliawiki (см. Границы надежности, основанные на достоверности матрицы Фишера). границы) и использовал описание для расчета дисперсии, а также верхней и нижней доверительных границ (R_U и R_L).

Вот пример рабочего кода для моей подгонки Вейбулла и моих доверительных границ с набором тестовых данных на основе описания 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. К сожалению, я действительно не понимаю, в чем проблема, и я не знаю, кого еще можно спросить. В инете и на разных форумах подходящего ответа не нашел.