From 58483c053d41186f3055b626fbe552f19cac3476 Mon Sep 17 00:00:00 2001 From: Matthew Nannemann Date: Sat, 20 Jun 2026 08:13:40 -0500 Subject: [PATCH 1/2] bug fix issue #57 --- reliability/Fitters.py | 434 ++++++++++++++++++++++------------------- tests/test_Fitters.py | 40 ++++ 2 files changed, 269 insertions(+), 205 deletions(-) diff --git a/reliability/Fitters.py b/reliability/Fitters.py index 4dc43956..d06c2c20 100644 --- a/reliability/Fitters.py +++ b/reliability/Fitters.py @@ -92,6 +92,105 @@ pd.options.display.width = 200 # prevents wrapping after default 80 characters +def _ensure_two_distinct_groups(group1, group2, model_name): + """ + Internal helper used by Fit_Weibull_Mixture and Fit_Weibull_CR. + + Final safety net: raises a clear ValueError (instead of a confusing inner + Weibull_2P error) if a group ends up with fewer than two distinct failure + times. The primary grouping logic (_get_group_dividing_line) now actively + tries to prevent this when the input data has sufficient distinct times. + """ + u1 = len(np.unique(group1)) + u2 = len(np.unique(group2)) + if u1 < 2 or u2 < 2: + raise ValueError( + f"The data provided does not appear to contain two distinct failure modes " + f"based on the automatic grouping algorithm used for the initial guess of the " + f"{model_name} model. This model could not be fitted. " + f"This commonly occurs with heavily right-censored data or when failures are " + f"heavily clustered. Consider the simpler distributions or exclude " + f"'{model_name}' when using Fit_Everything." + ) + + +def _get_group_dividing_line(failures): + """ + Internal helper used by Fit_Weibull_Mixture and Fit_Weibull_CR. + + Estimates the dividing line between the two groups for generating + initial parameter guesses. + + The base method fits a Gaussian KDE and draws chords from the peak + to the data range; the dividing line is chosen where the vertical + distance to the KDE is maximized. + + The result is then adjusted (when there are >=4 distinct failure times) + so that the split always produces at least two distinct failure times + in each group. This fixes the root cause of the old heuristic producing + degenerate subgroups even when a valid split existed among the unique + times. + """ + failures = np.asarray(failures) + if len(failures) == 0: + return 0.0 + + max_failures = np.max(failures) + min_failures = np.min(failures) + gkde = ss.gaussian_kde(failures) + delta = max_failures - min_failures + x_kde = np.linspace(min_failures - delta / 5, max_failures + delta / 5, 100) + y_kde = gkde.evaluate(x_kde) + peak_y = max(y_kde) + peak_x = x_kde[np.where(y_kde == peak_y)][0] + + left_x = min_failures + left_y = gkde.evaluate(left_x) + left_m = (peak_y - left_y) / (peak_x - left_x) + left_c = -left_m * left_x + left_y + left_line_x = np.linspace(left_x, peak_x, 100) + left_line_y = left_m * left_line_x + left_c # y=mx+c + left_kde = gkde.evaluate(left_line_x) + left_diff = abs(left_line_y - left_kde) + left_diff_max = max(left_diff) + left_div_line = left_line_x[np.where(left_diff == left_diff_max)][0] + + right_x = max_failures + right_y = gkde.evaluate(right_x) + right_m = (right_y - peak_y) / (right_x - peak_x) + right_c = -right_m * right_x + right_y + right_line_x = np.linspace(peak_x, right_x, 100) + right_line_y = right_m * right_line_x + right_c # y=mx+c + right_kde = gkde.evaluate(right_line_x) + right_diff = abs(right_line_y - right_kde) + right_diff_max = max(right_diff) + right_div_line = right_line_x[np.where(right_diff == right_diff_max)][0] + + if left_diff_max > right_diff_max: + dividing_line = left_div_line + else: + dividing_line = right_div_line + + # Adjust using distinct failure times (not just item counts) so that + # each side gets >= 2 distinct values whenever possible. + unique = np.unique(np.sort(failures)) + n_u = len(unique) + if n_u >= 4: + best_div = dividing_line + min_dist = float("inf") + # Possible splits that guarantee >=2 uniques on each side: + # split after index i (i >=1 and i <= n_u-3) + for i in range(1, n_u - 2): + cand = (unique[i] + unique[i + 1]) / 2.0 + dist = abs(cand - dividing_line) + if dist < min_dist: + min_dist = dist + best_div = cand + dividing_line = best_div + + return dividing_line + + class Fit_Everything: """ This function will fit all available distributions to the data provided. @@ -800,105 +899,119 @@ def __init__( ) if "Weibull_Mixture" not in self.excluded_distributions: - self.__Weibull_Mixture_params = Fit_Weibull_Mixture( - failures=failures, - right_censored=right_censored, - method=method, - optimizer=optimizer, - show_probability_plot=False, - print_results=False, - ) - self.Weibull_Mixture_alpha_1 = self.__Weibull_Mixture_params.alpha_1 - self.Weibull_Mixture_beta_1 = self.__Weibull_Mixture_params.beta_1 - self.Weibull_Mixture_alpha_2 = self.__Weibull_Mixture_params.alpha_2 - self.Weibull_Mixture_beta_2 = self.__Weibull_Mixture_params.beta_2 - self.Weibull_Mixture_proportion_1 = ( - self.__Weibull_Mixture_params.proportion_1 - ) - self.Weibull_Mixture_loglik = self.__Weibull_Mixture_params.loglik - self.Weibull_Mixture_BIC = self.__Weibull_Mixture_params.BIC - self.Weibull_Mixture_AICc = self.__Weibull_Mixture_params.AICc - self.Weibull_Mixture_AD = self.__Weibull_Mixture_params.AD - self.Weibull_Mixture_optimizer = self.__Weibull_Mixture_params.optimizer - self._parametric_CDF_Weibull_Mixture = ( - self.__Weibull_Mixture_params.distribution.CDF(xvals=d, show_plot=False) - ) - df = pd.concat( - [ - df, - pd.DataFrame( - data={ - "Distribution": ["Weibull_Mixture"], - "Alpha": [""], - "Beta": [""], - "Gamma": [""], - "Alpha 1": [self.Weibull_Mixture_alpha_1], - "Beta 1": [self.Weibull_Mixture_beta_1], - "Alpha 2": [self.Weibull_Mixture_alpha_2], - "Beta 2": [self.Weibull_Mixture_beta_2], - "Proportion 1": [self.Weibull_Mixture_proportion_1], - "DS": [""], - "Mu": [""], - "Sigma": [""], - "Lambda": [""], - "Log-likelihood": [self.Weibull_Mixture_loglik], - "AICc": [self.Weibull_Mixture_AICc], - "BIC": [self.Weibull_Mixture_BIC], - "AD": [self.Weibull_Mixture_AD], - "optimizer": [self.Weibull_Mixture_optimizer], - } - ), - ] - ) + try: + self.__Weibull_Mixture_params = Fit_Weibull_Mixture( + failures=failures, + right_censored=right_censored, + method=method, + optimizer=optimizer, + show_probability_plot=False, + print_results=False, + ) + self.Weibull_Mixture_alpha_1 = self.__Weibull_Mixture_params.alpha_1 + self.Weibull_Mixture_beta_1 = self.__Weibull_Mixture_params.beta_1 + self.Weibull_Mixture_alpha_2 = self.__Weibull_Mixture_params.alpha_2 + self.Weibull_Mixture_beta_2 = self.__Weibull_Mixture_params.beta_2 + self.Weibull_Mixture_proportion_1 = ( + self.__Weibull_Mixture_params.proportion_1 + ) + self.Weibull_Mixture_loglik = self.__Weibull_Mixture_params.loglik + self.Weibull_Mixture_BIC = self.__Weibull_Mixture_params.BIC + self.Weibull_Mixture_AICc = self.__Weibull_Mixture_params.AICc + self.Weibull_Mixture_AD = self.__Weibull_Mixture_params.AD + self.Weibull_Mixture_optimizer = self.__Weibull_Mixture_params.optimizer + self._parametric_CDF_Weibull_Mixture = ( + self.__Weibull_Mixture_params.distribution.CDF(xvals=d, show_plot=False) + ) + df = pd.concat( + [ + df, + pd.DataFrame( + data={ + "Distribution": ["Weibull_Mixture"], + "Alpha": [""], + "Beta": [""], + "Gamma": [""], + "Alpha 1": [self.Weibull_Mixture_alpha_1], + "Beta 1": [self.Weibull_Mixture_beta_1], + "Alpha 2": [self.Weibull_Mixture_alpha_2], + "Beta 2": [self.Weibull_Mixture_beta_2], + "Proportion 1": [self.Weibull_Mixture_proportion_1], + "DS": [""], + "Mu": [""], + "Sigma": [""], + "Lambda": [""], + "Log-likelihood": [self.Weibull_Mixture_loglik], + "AICc": [self.Weibull_Mixture_AICc], + "BIC": [self.Weibull_Mixture_BIC], + "AD": [self.Weibull_Mixture_AD], + "optimizer": [self.Weibull_Mixture_optimizer], + } + ), + ] + ) + except Exception as e: + colorprint( + f"WARNING: Weibull_Mixture was excluded from the results because it could not be fitted to the data.\n" + f"The error was: {str(e)[:200]}", + text_color="red", + ) if "Weibull_CR" not in self.excluded_distributions: - self.__Weibull_CR_params = Fit_Weibull_CR( - failures=failures, - right_censored=right_censored, - method=method, - optimizer=optimizer, - show_probability_plot=False, - print_results=False, - ) - self.Weibull_CR_alpha_1 = self.__Weibull_CR_params.alpha_1 - self.Weibull_CR_beta_1 = self.__Weibull_CR_params.beta_1 - self.Weibull_CR_alpha_2 = self.__Weibull_CR_params.alpha_2 - self.Weibull_CR_beta_2 = self.__Weibull_CR_params.beta_2 - self.Weibull_CR_loglik = self.__Weibull_CR_params.loglik - self.Weibull_CR_BIC = self.__Weibull_CR_params.BIC - self.Weibull_CR_AICc = self.__Weibull_CR_params.AICc - self.Weibull_CR_AD = self.__Weibull_CR_params.AD - self.Weibull_CR_optimizer = self.__Weibull_CR_params.optimizer - self._parametric_CDF_Weibull_CR = self.__Weibull_CR_params.distribution.CDF( - xvals=d, show_plot=False - ) - df = pd.concat( - [ - df, - pd.DataFrame( - data={ - "Distribution": ["Weibull_CR"], - "Alpha": [""], - "Beta": [""], - "Gamma": [""], - "Alpha 1": [self.Weibull_CR_alpha_1], - "Beta 1": [self.Weibull_CR_beta_1], - "Alpha 2": [self.Weibull_CR_alpha_2], - "Beta 2": [self.Weibull_CR_beta_2], - "Proportion 1": [""], - "DS": [""], - "Mu": [""], - "Sigma": [""], - "Lambda": [""], - "Log-likelihood": [self.Weibull_CR_loglik], - "AICc": [self.Weibull_CR_AICc], - "BIC": [self.Weibull_CR_BIC], - "AD": [self.Weibull_CR_AD], - "optimizer": [self.Weibull_CR_optimizer], - } - ), - ] - ) + try: + self.__Weibull_CR_params = Fit_Weibull_CR( + failures=failures, + right_censored=right_censored, + method=method, + optimizer=optimizer, + show_probability_plot=False, + print_results=False, + ) + self.Weibull_CR_alpha_1 = self.__Weibull_CR_params.alpha_1 + self.Weibull_CR_beta_1 = self.__Weibull_CR_params.beta_1 + self.Weibull_CR_alpha_2 = self.__Weibull_CR_params.alpha_2 + self.Weibull_CR_beta_2 = self.__Weibull_CR_params.beta_2 + self.Weibull_CR_loglik = self.__Weibull_CR_params.loglik + self.Weibull_CR_BIC = self.__Weibull_CR_params.BIC + self.Weibull_CR_AICc = self.__Weibull_CR_params.AICc + self.Weibull_CR_AD = self.__Weibull_CR_params.AD + self.Weibull_CR_optimizer = self.__Weibull_CR_params.optimizer + self._parametric_CDF_Weibull_CR = self.__Weibull_CR_params.distribution.CDF( + xvals=d, show_plot=False + ) + df = pd.concat( + [ + df, + pd.DataFrame( + data={ + "Distribution": ["Weibull_CR"], + "Alpha": [""], + "Beta": [""], + "Gamma": [""], + "Alpha 1": [self.Weibull_CR_alpha_1], + "Beta 1": [self.Weibull_CR_beta_1], + "Alpha 2": [self.Weibull_CR_alpha_2], + "Beta 2": [self.Weibull_CR_beta_2], + "Proportion 1": [""], + "DS": [""], + "Mu": [""], + "Sigma": [""], + "Lambda": [""], + "Log-likelihood": [self.Weibull_CR_loglik], + "AICc": [self.Weibull_CR_AICc], + "BIC": [self.Weibull_CR_BIC], + "AD": [self.Weibull_CR_AD], + "optimizer": [self.Weibull_CR_optimizer], + } + ), + ] + ) + except Exception as e: + colorprint( + f"WARNING: Weibull_CR was excluded from the results because it could not be fitted to the data.\n" + f"The error was: {str(e)[:200]}", + text_color="red", + ) if "Weibull_DS" not in self.excluded_distributions: self.__Weibull_DS_params = Fit_Weibull_DS( @@ -3611,10 +3724,15 @@ class Fit_Weibull_Mixture: than would be achieved by a Weibull mixture. For this reason, other types of mixtures are not implemented. - If the fitting process encounters a problem a warning will be printed. This - may be caused by the chosen distribution being a very poor fit to the data - or the data being heavily censored. If a warning is printed, consider trying - a different optimizer. + If the fitting process encounters a problem a warning will be printed (when + using Fit_Everything) or a clear ValueError will be raised (when calling the + fitter directly). Common causes include the data not containing two + distinguishable failure modes or the data being heavily censored. The + grouping heuristic used for the initial guess now chooses splits that + guarantee at least two distinct failures per subpopulation whenever the + overall data has four or more distinct failure times. A clear message is + produced that suggests excluding the model from Fit_Everything or trying the + simpler distributions. """ def __init__( @@ -3646,58 +3764,9 @@ def __init__( failures=failures, right_censored=right_censored ) # this is only used to find AD - # this algorithm is used to estimate the dividing line between the two groups - # firstly it fits a gaussian kde to the histogram - # then it draws two straight lines from the highest peak of the kde down to the lower and upper bounds of the failures - # the dividing line is the point where the difference between the kde and the straight lines is greatest - max_failures = max(failures) - min_failures = min(failures) - gkde = ss.gaussian_kde(failures) - delta = max_failures - min_failures - x_kde = np.linspace(min_failures - delta / 5, max_failures + delta / 5, 100) - y_kde = gkde.evaluate(x_kde) - peak_y = max(y_kde) - peak_x = x_kde[np.where(y_kde == peak_y)][0] - - left_x = min_failures - left_y = gkde.evaluate(left_x) - left_m = (peak_y - left_y) / (peak_x - left_x) - left_c = -left_m * left_x + left_y - left_line_x = np.linspace(left_x, peak_x, 100) - left_line_y = left_m * left_line_x + left_c # y=mx+c - left_kde = gkde.evaluate(left_line_x) - left_diff = abs(left_line_y - left_kde) - left_diff_max = max(left_diff) - left_div_line = left_line_x[np.where(left_diff == left_diff_max)][0] - - right_x = max_failures - right_y = gkde.evaluate(right_x) - right_m = (right_y - peak_y) / (right_x - peak_x) - right_c = -right_m * right_x + right_y - right_line_x = np.linspace(peak_x, right_x, 100) - right_line_y = right_m * right_line_x + right_c # y=mx+c - right_kde = gkde.evaluate(right_line_x) - right_diff = abs(right_line_y - right_kde) - right_diff_max = max(right_diff) - right_div_line = right_line_x[np.where(right_diff == right_diff_max)][0] - - if left_diff_max > right_diff_max: - dividing_line = left_div_line - else: - dividing_line = right_div_line - - number_of_items_in_group_1 = len(np.where(failures < dividing_line)[0]) - number_of_items_in_group_2 = len(failures) - number_of_items_in_group_1 - if number_of_items_in_group_1 < 2: - failures_sorted = np.sort(failures) - dividing_line = ( - failures_sorted[1] + failures_sorted[2] - ) / 2 # adjusts the dividing line in case there aren't enough failures in the first group - if number_of_items_in_group_2 < 2: - failures_sorted = np.sort(failures) - dividing_line = ( - failures_sorted[-2] + failures_sorted[-3] - ) / 2 # adjusts the dividing line in case there aren't enough failures in the second group + # Estimate dividing line (KDE heuristic + distinct-aware adjustment to + # guarantee >=2 distinct failures per group when the data supports it). + dividing_line = _get_group_dividing_line(failures) # this is the point at which data is assigned to one group or another for the purpose of generating the initial guess GROUP_1_failures = [] @@ -3715,6 +3784,8 @@ def __init__( else: GROUP_2_right_cens.append(item) + _ensure_two_distinct_groups(GROUP_1_failures, GROUP_2_failures, "Weibull_Mixture") + # 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, @@ -4144,58 +4215,9 @@ def __init__( failures=failures, right_censored=right_censored ) # this is only used to find AD - # this algorithm is used to estimate the dividing line between the two groups - # firstly it fits a gaussian kde to the histogram - # then it draws two straight lines from the highest peak of the kde down to the lower and upper bounds of the failures - # the dividing line is the point where the difference between the kde and the straight lines is greatest - max_failures = max(failures) - min_failures = min(failures) - gkde = ss.gaussian_kde(failures) - delta = max_failures - min_failures - x_kde = np.linspace(min_failures - delta / 5, max_failures + delta / 5, 100) - y_kde = gkde.evaluate(x_kde) - peak_y = max(y_kde) - peak_x = x_kde[np.where(y_kde == peak_y)][0] - - left_x = min_failures - left_y = gkde.evaluate(left_x) - left_m = (peak_y - left_y) / (peak_x - left_x) - left_c = -left_m * left_x + left_y - left_line_x = np.linspace(left_x, peak_x, 100) - left_line_y = left_m * left_line_x + left_c # y=mx+c - left_kde = gkde.evaluate(left_line_x) - left_diff = abs(left_line_y - left_kde) - left_diff_max = max(left_diff) - left_div_line = left_line_x[np.where(left_diff == left_diff_max)][0] - - right_x = max_failures - right_y = gkde.evaluate(right_x) - right_m = (right_y - peak_y) / (right_x - peak_x) - right_c = -right_m * right_x + right_y - right_line_x = np.linspace(peak_x, right_x, 100) - right_line_y = right_m * right_line_x + right_c # y=mx+c - right_kde = gkde.evaluate(right_line_x) - right_diff = abs(right_line_y - right_kde) - right_diff_max = max(right_diff) - right_div_line = right_line_x[np.where(right_diff == right_diff_max)][0] - - if left_diff_max > right_diff_max: - dividing_line = left_div_line - else: - dividing_line = right_div_line - - number_of_items_in_group_1 = len(np.where(failures < dividing_line)[0]) - number_of_items_in_group_2 = len(failures) - number_of_items_in_group_1 - if number_of_items_in_group_1 < 2: - failures_sorted = np.sort(failures) - dividing_line = ( - failures_sorted[1] + failures_sorted[2] - ) / 2 # adjusts the dividing line in case there aren't enough failures in the first group - if number_of_items_in_group_2 < 2: - failures_sorted = np.sort(failures) - dividing_line = ( - failures_sorted[-2] + failures_sorted[-3] - ) / 2 # adjusts the dividing line in case there aren't enough failures in the second group + # Estimate dividing line (KDE heuristic + distinct-aware adjustment to + # guarantee >=2 distinct failures per group when the data supports it). + dividing_line = _get_group_dividing_line(failures) # this is the point at which data is assigned to one group or another for the purpose of generating the initial guess GROUP_1_failures = [] @@ -4213,6 +4235,8 @@ def __init__( else: GROUP_2_right_cens.append(item) + _ensure_two_distinct_groups(GROUP_1_failures, GROUP_2_failures, "Weibull_CR") + # 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, diff --git a/tests/test_Fitters.py b/tests/test_Fitters.py index 96c375c5..00f65975 100644 --- a/tests/test_Fitters.py +++ b/tests/test_Fitters.py @@ -610,3 +610,43 @@ def test_Fit_Everything(): assert_allclose(LS.Exponential_1P_BIC, 197.06920107995282, rtol=rtol, atol=atol) assert_allclose(LS.Exponential_1P_loglik, -95.88544185670239, rtol=rtol, atol=atol) assert_allclose(LS.Exponential_1P_AD, 549.85986679373, rtol=rtol, atol=atol) + + +def test_Fit_Weibull_mixture_cr_robust_grouping(): + """Regression test: the grouping heuristic for Mixture/CR now uses distinct + failure times when adjusting the divider. Data with >=4 distinct failure + times (even one outlier + heavy clustering + heavy right-censoring at the + cluster) must produce valid (>=2 distinct per group) initial guesses and + not raise the "two distinct failure modes" error from the old heuristic. + """ + warnings.filterwarnings(action="ignore", category=RuntimeWarning) + + failed = [2] + [8] * 5 + [9] * 8 + [10] * 10 + passed = [10] * 75 + + # Direct calls should succeed (no distinct-group crash) + m = Fit_Weibull_Mixture( + failures=failed, right_censored=passed, + print_results=False, show_probability_plot=False + ) + assert m.alpha_1 > 0 and m.beta_1 > 0 + assert m.alpha_2 > 0 and m.beta_2 > 0 + assert m.proportion_1 > 0 + + c = Fit_Weibull_CR( + failures=failed, right_censored=passed, + print_results=False, show_probability_plot=False + ) + assert c.alpha_1 > 0 and c.beta_1 > 0 + assert c.alpha_2 > 0 and c.beta_2 > 0 + + # Fit_Everything should also succeed in fitting them (no exclusion) + fe = Fit_Everything( + failures=failed, right_censored=passed, + print_results=False, show_probability_plot=False, + show_histogram_plot=False, show_PP_plot=False, + show_best_distribution_probability_plot=False + ) + # Attributes exist (they were fitted) + assert fe.Weibull_Mixture_alpha_1 is not None + assert fe.Weibull_CR_alpha_1 is not None From 0ca7c5484cfc36fe15e08d44b0d585ad34243376 Mon Sep 17 00:00:00 2001 From: Matthew Nannemann Date: Sat, 20 Jun 2026 08:41:07 -0500 Subject: [PATCH 2/2] added changelog --- docs/Changelog.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/Changelog.rst b/docs/Changelog.rst index 3d0f5169..fb970c11 100644 --- a/docs/Changelog.rst +++ b/docs/Changelog.rst @@ -5,6 +5,14 @@ Changelog --------- +**Unreleased / next version** +''''''''''''''''''''''''''''' +**Bug Fixes** + +- Fixed a crash in ``Fit_Everything`` (and direct calls to ``Fit_Weibull_Mixture`` / ``Fit_Weibull_CR``) when the automatic two-group separation heuristic produced a subpopulation with fewer than two distinct failure times. Such data now produces a clear actionable ``ValueError`` (or is omitted from ``Fit_Everything`` results with a red warning). The original confusing inner "Weibull_2P needs 2 distinct failures" error is never seen by users. Addresses the repro case of one early outlier + heavy clustering at later times + heavy right censoring at the cluster value. A private helper ``_ensure_two_distinct_groups`` and try/except guards in ``Fit_Everything`` were added. Debug scaffolding left in the source during investigation was removed. + +- Root cause improvement to the grouping heuristic used by ``Fit_Weibull_Mixture`` and ``Fit_Weibull_CR`` for their initial guesses. The dividing line logic now operates on distinct failure times (instead of only raw item counts) and snaps to a split that guarantees at least two distinct failures per group whenever four or more distinct failure times are present in the data. This prevents the heuristic from unnecessarily producing degenerate subgroups. A shared helper ``_get_group_dividing_line`` was added (also reducing code duplication). The safety-net ``_ensure_two_distinct_groups`` remains in place. + **Version: 0.9.0 --- Released: 7 Mar 2025** '''''''''''''''''''''''''''''''''''''''''''