How to use the reliability.Fitters.Fit_Weibull_2P function in reliability

if type(failures) != np.ndarray:
            raise TypeError('failures must be a list or array of failure data')
        if type(right_censored) == list:
            right_censored = np.array(right_censored)
        if type(right_censored) != np.ndarray:
            raise TypeError('right_censored must be a list or array of right censored failure data')
        all_data = np.hstack([failures, right_censored])

        # solve it
        self.gamma = 0
        sp = ss.weibull_min.fit(all_data, floc=0, optimizer='powell')  # scipy's answer is used as an initial guess. Scipy is only correct when there is no censored data
        warnings.filterwarnings('ignore')  # necessary to supress the warning about the jacobian when using the nelder-mead optimizer

        if force_beta is None:
            guess = [sp[2], sp[0]]
            result = minimize(value_and_grad(Fit_Weibull_2P.LL), guess, args=(failures, right_censored), jac=True, tol=1e-6, method='nelder-mead')
        else:
            guess = [sp[2]]
            result = minimize(value_and_grad(Fit_Weibull_2P.LL_fb), guess, args=(failures, right_censored, force_beta), jac=True, tol=1e-6, method='nelder-mead')

        if result.success is True:
            params = result.x
            self.success = True
            if force_beta is None:
                self.alpha = params[0]
                self.beta = params[1]
            else:
                self.alpha = params * 1  # the *1 converts ndarray to float64
                self.beta = force_beta
        else:
            self.success = False
            print('WARNING: Fitting using Autograd FAILED for Weibull_2P. The fit from Scipy was used instead so results may not be accurate.')

right_censored = np.array(right_censored)
        if type(right_censored) != np.ndarray:
            raise TypeError('right_censored must be a list or array of right censored failure data')
        all_data = np.hstack([failures, right_censored])

        # solve it
        self.gamma = 0
        sp = ss.weibull_min.fit(all_data, floc=0, optimizer='powell')  # scipy's answer is used as an initial guess. Scipy is only correct when there is no censored data
        warnings.filterwarnings('ignore')  # necessary to supress the warning about the jacobian when using the nelder-mead optimizer

        if force_beta is None:
            guess = [sp[2], sp[0]]
            result = minimize(value_and_grad(Fit_Weibull_2P.LL), guess, args=(failures, right_censored), jac=True, tol=1e-6, method='nelder-mead')
        else:
            guess = [sp[2]]
            result = minimize(value_and_grad(Fit_Weibull_2P.LL_fb), guess, args=(failures, right_censored, force_beta), jac=True, tol=1e-6, method='nelder-mead')

        if result.success is True:
            params = result.x
            self.success = True
            if force_beta is None:
                self.alpha = params[0]
                self.beta = params[1]
            else:
                self.alpha = params * 1  # the *1 converts ndarray to float64
                self.beta = force_beta
        else:
            self.success = False
            print('WARNING: Fitting using Autograd FAILED for Weibull_2P. The fit from Scipy was used instead so results may not be accurate.')
            self.alpha = sp[2]
            self.beta = sp[0]

# within this loop, each list of failures and right censored values will be unpacked for each unique stress and plotted as a probability plot as well as the CDF of the common beta plot
        AICc_total = 0
        BIC_total = 0
        AICc = True
        for i, stress in enumerate(unique_stresses_f):
            failure_current_stress_df = f_df[f_df['stress'] == stress]
            FAILURES = failure_current_stress_df['times'].values
            if right_censored is not None:
                if stress in unique_stresses_rc:
                    right_cens_current_stress_df = rc_df[rc_df['stress'] == stress]
                    RIGHT_CENSORED = right_cens_current_stress_df['times'].values
                else:
                    RIGHT_CENSORED = None
            else:
                RIGHT_CENSORED = None
            weibull_fit_common_shape = Fit_Weibull_2P(failures=FAILURES, right_censored=RIGHT_CENSORED, show_probability_plot=False, print_results=False, force_beta=common_shape)
            weibull_fit_alpha_array_common_shape.append(weibull_fit_common_shape.alpha)
            if type(weibull_fit_common_shape.AICc) == str:
                AICc = False
            else:
                AICc_total += weibull_fit_common_shape.AICc
            BIC_total += weibull_fit_common_shape.BIC
            if show_plot is True:
                weibull_fit_common_shape.distribution.CDF(linestyle='--', color=color_list[i], xvals=xvals, plot_CI=False) # plotting of the confidence intervals has been turned off
                Probability_plotting.Weibull_probability_plot(failures=FAILURES, right_censored=RIGHT_CENSORED,plot_CI=False, color=color_list[i], label=str(stress))
                plt.legend(title='Stress')
                plt.xlim(10 ** (xmin + 1), 10 ** (xmax - 1))
                if common_shape_method == 'BIC':
                    plt.title(str('ALT Weibull Probability Plot\nOptimal BIC ' + r'$\beta$ = ' + str(round(common_shape, 4))))
                elif common_shape_method == 'weighted_average':
                    plt.title(str('ALT Weibull Probability Plot\nWeighted average ' + r'$\beta$ = ' + str(round(common_shape, 4))))
                elif common_shape_method == 'average':

def LL_fb(params, T_f, T_rc, force_beta):  # log likelihood function (2 parameter weibull) FORCED BETA
        LL_f = 0
        LL_rc = 0
        LL_f += Fit_Weibull_2P.logf(T_f, params[0], force_beta).sum()  # failure times
        LL_rc += Fit_Weibull_2P.logR(T_rc, params[0], force_beta).sum()  # right censored times
        return -(LL_f + LL_rc)

self.__Normal_2P_params = Fit_Normal_2P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
        self.Normal_2P_mu = self.__Normal_2P_params.mu
        self.Normal_2P_sigma = self.__Normal_2P_params.sigma
        self.Normal_2P_BIC = self.__Normal_2P_params.BIC
        self.Normal_2P_AICc = self.__Normal_2P_params.AICc
        self._parametric_CDF_Normal_2P = self.__Normal_2P_params.distribution.CDF(xvals=d, show_plot=False)

        self.__Lognormal_2P_params = Fit_Lognormal_2P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
        self.Lognormal_2P_mu = self.__Lognormal_2P_params.mu
        self.Lognormal_2P_sigma = self.__Lognormal_2P_params.sigma
        self.Lognormal_2P_BIC = self.__Lognormal_2P_params.BIC
        self.Lognormal_2P_AICc = self.__Lognormal_2P_params.AICc
        self._parametric_CDF_Lognormal_2P = self.__Lognormal_2P_params.distribution.CDF(xvals=d, show_plot=False)

        self.__Weibull_2P_params = Fit_Weibull_2P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
        self.Weibull_2P_alpha = self.__Weibull_2P_params.alpha
        self.Weibull_2P_beta = self.__Weibull_2P_params.beta
        self.Weibull_2P_BIC = self.__Weibull_2P_params.BIC
        self.Weibull_2P_AICc = self.__Weibull_2P_params.AICc
        self._parametric_CDF_Weibull_2P = self.__Weibull_2P_params.distribution.CDF(xvals=d, show_plot=False)

        self.__Gamma_2P_params = Fit_Gamma_2P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
        self.Gamma_2P_alpha = self.__Gamma_2P_params.alpha
        self.Gamma_2P_beta = self.__Gamma_2P_params.beta
        self.Gamma_2P_BIC = self.__Gamma_2P_params.BIC
        self.Gamma_2P_AICc = self.__Gamma_2P_params.AICc
        self._parametric_CDF_Gamma_2P = self.__Gamma_2P_params.distribution.CDF(xvals=d, show_plot=False)

        self.__Expon_1P_params = Fit_Expon_1P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
        self.Expon_1P_lambda = self.__Expon_1P_params.Lambda
        self.Expon_1P_BIC = self.__Expon_1P_params.BIC

def LL(params, T_f, T_rc):  # log likelihood function (2 parameter weibull)
        LL_f = 0
        LL_rc = 0
        LL_f += Fit_Weibull_2P.logf(T_f, params[0], params[1]).sum()  # failure times
        LL_rc += Fit_Weibull_2P.logR(T_rc, params[0], params[1]).sum()  # right censored times
        return -(LL_f + LL_rc)

params = [self.alpha, self.beta]
        k = len(params)
        n = len(all_data)
        LL2 = 2 * Fit_Weibull_2P.LL(params, failures, right_censored)
        self.loglik2 = LL2
        if n - k - 1 > 0:
            self.AICc = 2 * k + LL2 + (2 * k ** 2 + 2 * k) / (n - k - 1)
        else:
            self.AICc = 'Insufficient data'
        self.BIC = np.log(n) * k + LL2
        self.distribution = Weibull_Distribution(alpha=self.alpha, beta=self.beta)

        # confidence interval estimates of parameters
        Z = -ss.norm.ppf((1 - CI) / 2)
        if force_beta is None:
            hessian_matrix = hessian(Fit_Weibull_2P.LL)(np.array(tuple(params)), np.array(tuple(failures)), np.array(tuple(right_censored)))
            covariance_matrix = np.linalg.inv(hessian_matrix)
            self.alpha_SE = abs(covariance_matrix[0][0]) ** 0.5
            self.beta_SE = abs(covariance_matrix[1][1]) ** 0.5
            self.Cov_alpha_beta = abs(covariance_matrix[0][1])
            self.alpha_upper = self.alpha * (np.exp(Z * (self.alpha_SE / self.alpha)))
            self.alpha_lower = self.alpha * (np.exp(-Z * (self.alpha_SE / self.alpha)))
            self.beta_upper = self.beta * (np.exp(Z * (self.beta_SE / self.beta)))
            self.beta_lower = self.beta * (np.exp(-Z * (self.beta_SE / self.beta)))
        else:  # this is for when force beta is specified
            hessian_matrix = hessian(Fit_Weibull_2P.LL_fb)(np.array(tuple([self.alpha])), np.array(tuple(failures)), np.array(tuple(right_censored)), np.array(tuple([force_beta])))
            covariance_matrix = np.linalg.inv(hessian_matrix)
            self.alpha_SE = abs(covariance_matrix[0][0]) ** 0.5
            self.beta_SE = ''
            self.Cov_alpha_beta = ''
            self.alpha_upper = self.alpha * (np.exp(Z * (self.alpha_SE / self.alpha)))
            self.alpha_lower = self.alpha * (np.exp(-Z * (self.alpha_SE / self.alpha)))

xvals = np.logspace(-25, np.ceil(np.log10(max(failures))) + 6, 1000)

    if __fitted_dist_params is not None:
        if __fitted_dist_params.gamma > 0:
            fit_gamma = True

    if fit_gamma is False:
        if __fitted_dist_params is not None:
            alpha = __fitted_dist_params.alpha
            beta = __fitted_dist_params.beta
            alpha_SE = __fitted_dist_params.alpha_SE
            beta_SE = __fitted_dist_params.beta_SE
            Cov_alpha_beta = __fitted_dist_params.Cov_alpha_beta
        else:
            from reliability.Fitters import Fit_Weibull_2P
            fit = Fit_Weibull_2P(failures=failures, right_censored=right_censored, CI=CI, show_probability_plot=False, print_results=False)
            alpha = fit.alpha
            beta = fit.beta
            alpha_SE = fit.alpha_SE
            beta_SE = fit.beta_SE
            Cov_alpha_beta = fit.Cov_alpha_beta
        if 'label' in kwargs:
            label = kwargs.pop('label')
        else:
            label = str('Fitted Weibull_2P (α=' + str(round_to_decimals(alpha, dec)) + ', β=' + str(round_to_decimals(beta, dec)) + ')')
        if 'color' in kwargs:
            data_color = kwargs.get('color')
        else:
            data_color = 'k'
        xlabel = 'Time'
    elif fit_gamma is True:
        if __fitted_dist_params is not None:

'''
            __BIC_minimizer is used by the minimize function to get the shape which gives the lowest overall BIC
            '''
            BIC_tot = 0
            for stress in unique_stresses_f:
                failure_current_stress_df = f_df[f_df['stress'] == stress]
                FAILURES = failure_current_stress_df['times'].values
                if right_censored is not None:
                    if stress in unique_stresses_rc:
                        right_cens_current_stress_df = rc_df[rc_df['stress'] == stress]
                        RIGHT_CENSORED = right_cens_current_stress_df['times'].values
                    else:
                        RIGHT_CENSORED = None
                else:
                    RIGHT_CENSORED = None
                weibull_fit_common_shape = Fit_Weibull_2P(failures=FAILURES, right_censored=RIGHT_CENSORED, show_probability_plot=False, print_results=False, force_beta=common_shape_X)
                BIC_tot += weibull_fit_common_shape.BIC
            return BIC_tot

GROUP_2_failures = []
        GROUP_1_right_cens = []
        GROUP_2_right_cens = []
        for item in failures:
            if item < division_line:
                GROUP_1_failures.append(item)
            else:
                GROUP_2_failures.append(item)
        for item in right_censored:
            if item < division_line:
                GROUP_1_right_cens.append(item)
            else:
                GROUP_2_right_cens.append(item)

        # get inputs for the guess by fitting a weibull to each of the groups with their respective censored data
        group_1_estimates = Fit_Weibull_2P(failures=GROUP_1_failures, right_censored=GROUP_1_right_cens, show_probability_plot=False, print_results=False)
        group_2_estimates = Fit_Weibull_2P(failures=GROUP_2_failures, right_censored=GROUP_2_right_cens, show_probability_plot=False, print_results=False)
        p_guess = (len(GROUP_1_failures) + len(GROUP_1_right_cens)) / len(all_data)  # proportion guess
        guess = [group_1_estimates.alpha, group_1_estimates.beta, group_2_estimates.alpha, group_2_estimates.beta, p_guess]  # A1,B1,A2,B2,P

        # solve it
        bnds = [(0.0001, None), (0.0001, None), (0.0001, None), (0.0001, None), (0.0001, 0.9999)]  # bounds of solution
        result = minimize(value_and_grad(Fit_Weibull_Mixture.LL), guess, args=(failures, right_censored), jac=True, bounds=bnds, tol=1e-6)
        params = result.x
        self.alpha_1 = params[0]
        self.beta_1 = params[1]
        self.alpha_2 = params[2]
        self.beta_2 = params[3]
        self.proportion_1 = params[4]
        self.proportion_2 = 1 - params[4]

        params = [self.alpha_1, self.beta_1, self.alpha_2, self.beta_2, self.proportion_1]