.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/reliability/life_data/plot_02_cohen_censored_mle.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_reliability_life_data_plot_02_cohen_censored_mle.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_reliability_life_data_plot_02_cohen_censored_mle.py:


Cohen's maximum likelihood for censored Weibull data (1965)
===========================================================

A rank-regression line uses a censored unit only through the ranks of
later failures. A. Clifford Cohen (1965) wrote out the maximum-likelihood
equations of the Weibull distribution for complete, singly censored and
progressively censored samples. Every unit contributes: a failure through
its density, a survivor through :math:`R(t)`. The shape solves one
equation in :math:`\beta` alone, and :math:`\eta` follows
(:func:`~engineeringkit.reliability.weibull_mle`). The example repeats a
heavily censored test many times and compares the scatter of the two
estimators.

.. GENERATED FROM PYTHON SOURCE LINES 17-22

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from engineeringkit.reliability import weibull_mle, weibull_rank_regression








.. GENERATED FROM PYTHON SOURCE LINES 23-28

One test, stopped at 600 h
--------------------------

Twenty units with :math:`\eta` = 1000 h, :math:`\beta` = 1.8; the test
ends at 600 h with the survivors still running (Type I censoring).

.. GENERATED FROM PYTHON SOURCE LINES 28-44

.. code-block:: Python


    eta, beta, n, stop = 1000.0, 1.8, 20, 600.0
    rng = np.random.default_rng(7)


    def run_test():
        life = eta * rng.weibull(beta, n)
        return np.minimum(life, stop), (life <= stop).astype(int)


    times, failed = run_test()
    mle, rr = weibull_mle(times, failed), weibull_rank_regression(times, failed)
    print(f"{failed.sum()} of {n} failed by {stop:.0f} h")
    print(f"maximum likelihood: beta = {mle.beta:.2f}, eta = {mle.eta:.0f} h")
    print(f"rank regression:    beta = {rr.beta:.2f}, eta = {rr.eta:.0f} h")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    6 of 20 failed by 600 h
    maximum likelihood: beta = 2.07, eta = 1001 h
    rank regression:    beta = 0.97, eta = 2742 h




.. GENERATED FROM PYTHON SOURCE LINES 45-47

Scatter over 2000 repeated tests
--------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 47-71

.. code-block:: Python


    estimates = {"maximum likelihood": [], "rank regression": []}
    while len(estimates["rank regression"]) < 2000:
        times, failed = run_test()
        if failed.sum() >= 3:
            estimates["maximum likelihood"].append(weibull_mle(times, failed).beta)
            estimates["rank regression"].append(weibull_rank_regression(times, failed).beta)
    fig, ax = plt.subplots(figsize=(6, 3.8))
    for label, values in estimates.items():
        values = np.array(values)
        rmse = np.sqrt(np.mean((values - beta) ** 2))
        print(f"{label:18s}: median beta {np.median(values):.2f}, RMS error {rmse:.2f}")
        ax.hist(values, bins=np.linspace(0, 5, 60), histtype="step", label=f"{label} (RMS error {rmse:.2f})")
    ax.axvline(beta, color="k", ls=":")
    ax.set(xlabel=r"estimated $\beta$", ylabel="tests", title="20 units, test stopped at 600 h")
    ax.legend()
    fig.tight_layout()

    mle_error = np.sqrt(np.mean((np.array(estimates["maximum likelihood"]) - beta) ** 2))
    rr_error = np.sqrt(np.mean((np.array(estimates["rank regression"]) - beta) ** 2))
    assert abs(np.median(estimates["maximum likelihood"]) - beta) < 0.2
    assert mle_error < rr_error

    plt.show()



.. image-sg:: /api/gallery/reliability/life_data/images/sphx_glr_plot_02_cohen_censored_mle_001.png
   :alt: 20 units, test stopped at 600 h
   :srcset: /api/gallery/reliability/life_data/images/sphx_glr_plot_02_cohen_censored_mle_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    maximum likelihood: median beta 1.91, RMS error 1.13
    rank regression   : median beta 1.70, RMS error 2.04





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.729 seconds)


.. _sphx_glr_download_api_gallery_reliability_life_data_plot_02_cohen_censored_mle.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../../lite/lab/index.html?path=api/gallery/reliability/life_data/plot_02_cohen_censored_mle.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_02_cohen_censored_mle.ipynb <plot_02_cohen_censored_mle.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_02_cohen_censored_mle.py <plot_02_cohen_censored_mle.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_02_cohen_censored_mle.zip <plot_02_cohen_censored_mle.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_