.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/seismology/plot_11_gutenberg_richter_law.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_seismology_plot_11_gutenberg_richter_law.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_seismology_plot_11_gutenberg_richter_law.py:


Gutenberg and Richter's frequency-magnitude law (1944)
======================================================

Beno Gutenberg and Charles Richter found that the number of earthquakes of at
least magnitude M falls as 10^(a - bM), with b close to 1: each magnitude unit
brings ten times fewer events. A synthetic catalogue drawn from that law
recovers b with Aki's maximum-likelihood estimator.

.. GENERATED FROM PYTHON SOURCE LINES 10-33



.. image-sg:: /api/gallery/seismology/images/sphx_glr_plot_11_gutenberg_richter_law_001.png
   :alt: Estimated b = 1.003
   :srcset: /api/gallery/seismology/images/sphx_glr_plot_11_gutenberg_richter_law_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from earthsciencekit.seismology import gutenberg_richter

    rng = np.random.default_rng(1944)
    b, m_min, n = 1.0, 2.0, 20000
    magnitudes = m_min + rng.exponential(1 / (b * np.log(10)), n)
    a = np.log10(n) + b * m_min
    b_hat = np.log10(np.e) / (magnitudes.mean() - m_min)
    assert abs(b_hat - b) < 0.03
    assert np.isclose(gutenberg_richter(a, b, m_min), n)

    m = np.linspace(m_min, magnitudes.max(), 200)
    observed = (magnitudes[None, :] >= m[:, None]).sum(axis=1)
    fig, ax = plt.subplots()
    ax.semilogy(m, observed, ".", ms=3, label="synthetic catalogue")
    ax.semilogy(m, gutenberg_richter(a, b, m), label=f"law with b = {b}")
    ax.set(xlabel="magnitude M", ylabel="number with magnitude >= M", title=f"Estimated b = {b_hat:.3f}")
    ax.legend()

    plt.show()


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

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


.. _sphx_glr_download_api_gallery_seismology_plot_11_gutenberg_richter_law.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/seismology/plot_11_gutenberg_richter_law.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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