updating reliability_growth and CI docs

docs/Changelog.rst

@@ -10,13 +10,16 @@ Changelog
 
 **Summary of changes**
 
-This will be writted closer to the release date.
+This will be written closer to the release date.
 
 **New features**
 
+-    New dataset called system_growth has been added to the Datasets module.
+-    Repairable_systems.reliability growth has been completely rewritten. This function now includes both the Duane and Crow-AMSAA reliability growth models. The parametrisation of the Duane model has been modified to match what reliasoft uses.
 
 **API Changes**
 
+-    Repairable_systems.reliability growth has been completely rewritten, so the inputs and outputs are completely different. The old version of the Duane model has been replaced without deprecation. Users needing to use the old version should use v0.7.1 of reliability.
 
 **Bug Fixes**
 
@@ -26,6 +29,8 @@ This will be writted closer to the release date.
 
 -    Improvements to the API documentation for Convert_data, Datasets, PoF, and Utils modules. This has been a long term body of work to reformat the documentation, and it is finally complete.
 -    The required version of matplotlib has been upgraded to 3.5.0 to enable the above bugfix for the computed_zorder in ALT life stress plots.
+-    Theory documents are finished for `censored data <https://reliability.readthedocs.io/en/latest/What%20is%20censored%20data.html>`_, `plotting positions <https://reliability.readthedocs.io/en/latest/How%20are%20the%20plotting%20positions%20calculated.html>`_, `Least Squares Estimation <https://reliability.readthedocs.io/en/latest/How%20does%20Least%20Squares%20Estimation%20work.html>`_, `Maximum Likelihood Estimation <https://reliability.readthedocs.io/en/latest/How%20does%20Maximum%20Likelihood%20Estimation%20work.html>`_, and `Confidence Intervals <https://reliability.readthedocs.io/en/latest/How%20are%20the%20confidence%20intervals%20calculated.html>`_.
+-    Updates pytests for new reliability_growth function.
 
 **Version: 0.7.1 --- Released: 26 Oct 2021**
 ''''''''''''''''''''''''''''''''''''''''''''

docs/Datasets.rst

@@ -15,6 +15,7 @@ There are a few datasets that have been included with reliability that users may
 
 - automotive - 10 failures, 21 right censored. It is used in `this example <https://reliability.readthedocs.io/en/latest/Kaplan-Meier%20estimate%20of%20reliability.html>`_
 - mileage - 100 failures with no right censoring. It is used in the examples for `KStest <https://reliability.readthedocs.io/en/latest/Kolmogorov-Smirnov%20test.html>`_ and `chi2test <https://reliability.readthedocs.io/en/latest/Chi-squared%20test.html>`_.
+- system_growth - 22 failures with no right censoring. It is used in the example for `reliability_growth <https://reliability.readthedocs.io/en/latest/Reliability%20growth.html>`_.
 - defective_sample - 1350 failures, 12296 right censored. It exhibits the behavior of a defective sample (also known as Limited failure population or Defective subpopulation).
 - mixture - 71 failures, 3320 right censored. This is best modelled using a mixture model.
 - electronics - 10 failures, 4072 right censored. It is used in `this example <https://reliability.readthedocs.io/en/latest/Fitting%20a%20specific%20distribution%20to%20data.html#using-fit-weibull-2p-grouped-for-large-data-sets>`_.

docs/How are the confidence intervals calculated.rst

@@ -5,9 +5,9 @@
 How are the confidence intervals calculated
 '''''''''''''''''''''''''''''''''''''''''''
 
-There confidence intervals on the model parameters, and confidence intervals on the plots.
+There are confidence intervals on the model parameters, and confidence intervals on the plots. This document explains how both are calculated.
 The confidence intervals on the plots use the confidence intervals on the model parameters in the procedure that is explained in the second section of this document.
-Before reading this document, you should have a reasonable understanding of how `Maximum Likelihood Estimation (MLE) works <https://reliability.readthedocs.io/en/latest/How%20does%20Maximum%20Likelihood%20Estimation%20work.html>`_ and also understand what `partial derivatives <https://www.derivative-calculator.net/>`_ are.
+Before reading this document, you should have a reasonable understanding of how `Maximum Likelihood Estimation (MLE) works <https://reliability.readthedocs.io/en/latest/How%20does%20Maximum%20Likelihood%20Estimation%20work.html>`_ and also understand `what partial derivatives are <https://en.wikipedia.org/wiki/Partial_derivative>`_ and `how to calculate them <https://www.derivative-calculator.net/>`_.
 
 Confidence intervals on the parameters
 """"""""""""""""""""""""""""""""""""""
@@ -243,13 +243,13 @@ In code (for bounds on time) it looks like this:
             + 2 * dv_da(u, self.alpha, self.beta) * dv_db(u, self.alpha, self.beta) * self.Cov_alpha_beta
         )
 
-     # v is ln(t) and Y is reliability
-     v_lower = v(Y, self.alpha, self.beta) - Z * (var_v(self, Y) ** 0.5)
-     v_upper = v(Y, self.alpha, self.beta) + Z * (var_v(self, Y) ** 0.5)
+    # v is ln(t) and Y is reliability
+    v_lower = v(Y, self.alpha, self.beta) - Z * (var_v(self, Y) ** 0.5)
+    v_upper = v(Y, self.alpha, self.beta) + Z * (var_v(self, Y) ** 0.5)
 
-     # transform back from v = ln(t)
-     t_lower = np.exp(v_lower)
-     t_upper = np.exp(v_upper)
+    # transform back from v = ln(t)
+    t_lower = np.exp(v_lower)
+    t_upper = np.exp(v_upper)
 
 There are several other challenges to getting Python to do this correctly, such as where to incorporate :math:`\gamma` for location shifted distributions, how to distribute the points so they look smooth, how to correct for things (like reversals in the bounds) that are mathematically correct but practically (in the world of reliability engineering) incorrect, and how to correctly transform the bounds on the SF to get the bounds on the CDF or CHF.
 Some distributions (such as the Gamma Distribution) are particularly difficult and require a slightly different method to that which is explained above.
@@ -281,7 +281,7 @@ For bounds on time, this can be done programatically using the percentiles optio
 
 .. image:: images/CI_example2.png
 
-For bounds on reliability, extracting the parameters programatically is not currently enabled. It will be part of a future release (likely in January 2021), and will be available directly from the plotting method (avoiding the complicated ways shown above).
+For bounds on reliability, extracting the parameters programatically is not currently enabled. It will be part of a future release (likely in January 2021), and will be available directly from the plotting method (avoiding the complicated method shown above).
 
 For bounds on either time or reliability, the ``Other_functions.crosshairs`` function provides an interactive set of crosshairs which can be used to find the values using the mouse.
 A demo of how this works is shown `here <https://reliability.readthedocs.io/en/latest/Crosshairs.html>`_.

docs/Reliability growth.rst

@@ -5,30 +5,180 @@
 Reliability growth
 ''''''''''''''''''
 
-The reliability of a non-repairable component always decreases with time, but for repairable systems the term "reliability growth" refers to the process of gradual product improvement through the elimination of design deficiencies. In repairable systems, reliability growth is observable through an increase in the interarrival times of failures. Reliability growth is applicable to all levels of design decomposition from complete systems down to components. The maximum achieveable reliability is locked in by design, so reliability growth above the design reliability is only possible through design changes. It may be possible to achieve some reliability growth through other improvements (such as to the maintenance program) though these improvements will only help the system to achieve its design reliability.
+The reliability of a non-repairable component always decreases with time, but for repairable systems the term "reliability growth" refers to the process of gradual product improvement through the elimination of design deficiencies.
+In repairable systems, reliability growth is observable through an increase in the interarrival times of failures.
+Reliability growth is applicable to all levels of design decomposition from complete systems down to components.
+The maximum achieveable reliability is locked in by design, so reliability growth above the design reliability is only possible through design changes.
+It may be possible to achieve some reliability growth through other improvements (such as optimizing the maintenance program) though these improvements will only help the system to achieve its design reliability.
 
-The Duane method of modeling reliability growth involves the use of the total time on test [t] (we may also use distance, cycles, etc.) when the failure occurred and the sequence of the failure [N]. The cumulative mean time between failures (MTBF) is :math:`t/N`. By plotting :math:`ln(t)` vs :math:`ln(t/N)` we obtain a straight line which is used get the parameters lambda and beta. Using these parameters, we can model the instantaneous MTBF in the form :math:`\frac{t^{1-\beta}}{\lambda \times \beta}`. The function `reliability_growth` accepts the failure times and performs this model fitting to obtain the parameters lambda and beta, as well as produce the reliability growth plot. It is often of interest to know how much total time on test we need to meet a target MTBF. This can be found analytically by specifying the target_MTBF argument.
+The function `Repairable_systems.reliability_growth` fits a model to the failure times so that the growth (or deterioration) of the mean time between failures (MTBF) can be predicted.
+The two models available within `reliability` are the Duane model and the Crow-AMSAA model.
+
+In both the Duane and Crow-AMSAA models, the input is a list of failure times (actual times not interarrival times).
+The point estimate for the cumulative MTBF is calculated as :math:`\frac{{\rm failure time}}{{\rm failure number}}`.
+This produces the scatter plot shown in the plots below.
+
+The smooth curve shows the model (Duane or Crow-AMSAA) that is fitted to the data. The forumlation of these models is explained below.
+
+Duane model
+"""""""""""
+
+The algorithm to fit the model is as follows:
+
+1. accept an input of the failures times and sort the times in ascending order. Let the largest time be :math:`T`.
+2. create an array of failure numbers from 1 to n.
+3. calculate :math:`MTBF_{cumulative} = \frac{{\rm failure times}}{{\rm failure number}}`. This is the scatter points seen on the plot.
+4. Convert to log space: :math:`x = ln({\rm failure times})` and :math:`y = ln(MTBF_{cumulative})`
+5. fit a straight line (:math:`y=m.x+c`) to the data to get the model parameters.
+6. extract the model parameters from the parameters of the straight line, such that :math:`\alpha = m` and :math:`b = exp(c)`
+
+This gives us the model parameters of :math:`b` and :math:`\alpha`. The formulas for the other reported values are:
+
+:math:`DMTBF_C = b \times (T^{\alpha})`. This is the demonstrated MTBF (cumulative) and is reported in the results as DMTBF_C.
+
+:math:`DFI_C = \frac{1}{DMTBF_C}`. This is the demonstrated failure intensity (cumulative) and is reported in the results as DFI_C.
+
+:math:`DFI_I = (1 - \alpha) \times DFI_C`. This is the demonstrated failure intensity (instantaneous) and is reported in the results as DFI_I.
+
+:math:`DMTBF_I = \frac{1}{DFI_I}`. This is the demonstrated MTBF (instantaneous) and is reported in the results as DMTBF_I.
+
+:math:`A = \frac{1}{b}`. This is reported in the results as A.
+
+The time to reach the target MTBF is calculated as :math:`t_{target} = \left( \frac{target_{MTBF}}{b} \right)^{\frac{1}{\alpha}}`
+
+For more information see `reliawiki <http://reliawiki.org/index.php/Duane_Model>`_.
+
+Crow-AMSAA model
+""""""""""""""""
+
+The algorithm to fit the model is as follows:
+
+1. accept an input of the failures times and sort the times in ascending order. Let the largest time be :math:`T`.
+2. create an array of failure numbers from 1 to n.
+3. calculate :math:`MTBF_{cumulative} = \frac{{\rm failure times}}{{\rm failure number}}`. This is the scatter points seen on the plot.
+4. calculate :math:`\beta = \frac{n}{n \times ln(T) - \sum ln({\rm failure times})}`
+5. calculate :math:`\lambda = \frac{n}{T^{\beta}}`
+
+This gives us the model parameters :math:`\beta` and :math:`\lambda`. The formulas for the other reported values are:
+
+:math:`growth rate = 1 - \beta`. This is reported in the results as growth_rate.
+
+:math:`DMTBF_I = \frac{1}{\lambda \times \beta \times T^{\beta - 1}}`. This is the demonstrated MTBF (instantaneous) and is reported in the results as DMTBF_I.
+
+:math:`DFI_I = \frac{1}{DMTBF_I}`. This is the demonstrated failure intensity (instantaneous) and is reported in the results as DFI_I.
+
+:math:`DMTBF_C = \frac{1}{\lambda} \times T^{1 - \beta}`. This is the demonstrated MTBF (cumulative) and is reported in the results as DMTBF_C.
+
+:math:`DFI_C = \frac{1}{DMTBF_C}`. This is the demonstrated failure intensity (cumulative) and is reported in the results as DFI_C.
+
+The time to reach the target MTBF is calculated as :math:`t_{target} = \left(\frac{1}{\lambda \times \target_{MTBF}} \right)^ \frac{1}{\beta - 1}`
+
+For more information see `reliawiki <http://reliawiki.org/index.php/Crow-AMSAA_(NHPP)>`_.
 
 .. admonition:: API Reference
 
    For inputs and outputs see the `API reference <https://reliability.readthedocs.io/en/latest/API/Repairable_systems/reliability_growth.html>`_.
 
-In the example below, we supply the total time on test when each failure occurred, and we also supply the target_MTBF as 50000 so that we can find out how much total time on test will be needed to reach the target MTBF.
+Example 1
+"""""""""
+
+In this first example, we import a dataset and fit the Duane model. Log_scale is set to True (which only affects the plot). The target MTBF is 35 which will give us the time to reach the target MTBF based on the model.
 
 .. code:: python
 
     from reliability.Repairable_systems import reliability_growth
+    from reliability.Datasets import system_growth
     import matplotlib.pyplot as plt
-    times = [10400,26900,43400,66400,89400,130400,163400,232000,242000,340700]
-    reliability_growth(times=times,target_MTBF=50000,label='Reliability growth curve',xmax=500000)
-    plt.legend()
+    reliability_growth(times=system_growth().failures, model="Duane", target_MTBF=35, log_scale=True)
     plt.show()
+
+    '''
+    Duane Reliability growth model parameters:
+    Alpha: 0.42531
+    A: 0.57338
+    Demonstrated MTBF (cumulative): 26.86511
+    Demonstrated MTBF (instantaneous): 46.74719
+    Demonstrated failure intensity (cumulative): 0.037223
+    Demonstrated failure intensity (instantaneous): 0.021392
+    Time to reach target MTBF: 1154.77862
+    '''
+
+.. image:: images/reliability_growth_Duane.png
+
+Example 2
+"""""""""
+
+In this second example, we are using the same failure times as the example above, but now we are fitting the Crow-AMSAA model. The MTBF plot is in linear scale since log_scale has not been specified and it defaults to False.
+Once again, the target MTBF of 35 is specified and the results tell us the time to reach this target.
+
+.. code:: python
+
+    from reliability.Repairable_systems import reliability_growth
+    from reliability.Datasets import system_growth
+    import matplotlib.pyplot as plt
+
+    reliability_growth(times=system_growth().failures, model="Crow-AMSAA", target_MTBF=35)
+    plt.show()
+
+    '''
+    Crow-AMSAA Reliability growth model parameters:
+    Beta: 0.61421
+    Lambda: 0.42394
+    Growth rate: 0.38579
+    Demonstrated MTBF (cumulative): 28.18182
+    Demonstrated MTBF (instantaneous): 45.883
+    Demonstrated failure intensity (cumulative): 0.035484
+    Demonstrated failure intensity (instantaneous): 0.021795
+    Time to reach target MTBF: 1087.18769
+    '''
+
+.. image:: images/reliability_growth_CrowAMSAA.png
+
+Example 3
+"""""""""
+
+In this third example, we will compare the two models in both linear space (left plot) and log space (right plot).
+The fit of the Duane model through the points seems much better than is achieved by the Crow-AMSAA model, though this depends on the dataset.
+The Crow-AMSAA model places a strong emphasis on the last data point and will always ensure the model passes through this point.
+Depending on whether the last data point sits above or below the average will affect whether the Crow-AMSAA model is more optimistic (higher) or pessimistic (lower) in its prediction of the achieved MTBF than that which is predicted by the Duane model.
+
+.. code:: python
+
+    from reliability.Repairable_systems import reliability_growth
+    from reliability.Datasets import automotive
+    import matplotlib.pyplot as plt
     
+    plt.figure(figsize=(10,5))
+    plt.subplot(121)
+    reliability_growth(times=automotive().failures, model="Crow-AMSAA", target_MTBF=14000)
+    reliability_growth(times=automotive().failures, model="Duane", target_MTBF=14000,color='darkorange')
+    plt.subplot(122)
+    reliability_growth(times=automotive().failures, model="Crow-AMSAA", target_MTBF=14000,print_results=False,log_scale=True)
+    reliability_growth(times=automotive().failures, model="Duane", target_MTBF=14000,color='darkorange',print_results=False,log_scale=True)
+    plt.title('MTBF vs Time (log scale)')
+    plt.show()
+
     '''
-    Reliability growth model parameters:
-    lambda: 0.002355878294089656
-    beta: 0.6638280053477188
-    Time to reach target MTBF: 428131.18039448344
+    Crow-AMSAA Reliability growth model parameters:
+    Beta: 0.67922
+    Lambda: 0.0033282
+    Growth rate: 0.32078
+    Demonstrated MTBF (cumulative): 13190
+    Demonstrated MTBF (instantaneous): 19419.22019
+    Demonstrated failure intensity (cumulative): 7.5815e-05
+    Demonstrated failure intensity (instantaneous): 5.1495e-05
+    Time to reach target MTBF: 158830.62457
+    
+    Duane Reliability growth model parameters:
+    Alpha: 0.3148
+    A: 0.0038522
+    Demonstrated MTBF (cumulative): 10620.71841
+    Demonstrated MTBF (instantaneous): 15500.20608
+    Demonstrated failure intensity (cumulative): 9.4156e-05
+    Demonstrated failure intensity (instantaneous): 6.4515e-05
+    Time to reach target MTBF: 317216.14347
     '''
 
-.. image:: images/reliability_growth.png
+.. image:: images/reliability_growth_both.png
+
+.. note:: The function ``reliability_growth`` was completely rewritten in v0.8.0 to match the method used by Reliasoft. Prior to v0.8.0, only the Duane model was available, and the values returned were for a model with a completely different parameterisation.

reliability/Datasets.py

@@ -114,6 +114,50 @@ def __init__(self):
         self.info = df
 
 
+class system_growth:
+    """
+    This dataset is contains 22 values with no right censoring.
+    The data is from a system that has an increasing MTBF.
+
+    Returns
+    -------
+    failures : list
+        The failure times
+    info : dataframe
+        Descriptive statistics about the dataset
+
+    Notes
+    -----
+    When importing the dataset, ensure it is called using the brackets after the
+    name. Example Usage:
+
+    .. code:: python
+
+        from reliability.Datasets import system_growth
+        from reliability.Fitters import Fit_Weibull_2P
+        Fit_Weibull_2P(failures=system_growth().failures)
+    """
+
+    def __init__(self):
+        self.failures = [2.7, 10.3, 12.5, 30.6, 57, 61.3, 80, 109.5, 125, 128.6, 143.8, 167.9, 229.2, 296.7, 320.6, 328.2, 366.2, 396.7, 421.1, 438.2, 501.2, 620]
+        rc = 0
+        f = len(self.failures)
+        tot = f + rc
+        data = {
+            "Stat": ["Name", "Total Values", "Failures", "Right Censored"],
+            "Value": [
+                "system_growth",
+                tot,
+                str(str(f) + " (" + str(round(f / tot * 100, 2)) + "%)"),
+                str(str(rc) + " (" + str(round(rc / tot * 100, 2)) + "%)"),
+            ],
+        }
+        df = pd.DataFrame(data, columns=["Stat", "Value"])
+        blankIndex = [""] * len(df)
+        df.index = blankIndex
+        self.info = df
+
+
 class defective_sample:
     """
     This dataset is heavily right censored with intermixed multiply censored