Merge branch 'main' of github.com:swiss-seismological-service/SeismoS…

…tats into tests/catalog-methods

docs/source/reference/analysis.md

@@ -4,6 +4,15 @@
 .. currentmodule:: seismostats
 ```
 
+## Estimating the Completeness Magnitude
+```{eval-rst}
+.. autosummary::
+    :toctree: api/
+
+    analysis.mc_ks
+    analysis.mc_max_curvature
+```
+
 ## Estimating b-Values
 
 ```{eval-rst}
@@ -12,7 +21,7 @@
 
     estimate_b
     shi_bolt_confidence
-    analysis.estimate_b_tinti
+    analysis.estimate_b_classic
     analysis.estimate_b_positive
     analysis.estimate_b_more_positive
     analysis.estimate_b_utsu
@@ -21,13 +30,15 @@
     analysis.estimate_b_kijko_smit
 ```
 
-## Estimating the Completeness Magnitude
+## Estimating a-Values
+
 ```{eval-rst}
 .. autosummary::
     :toctree: api/
 
-    analysis.mc_ks
-    analysis.mc_max_curvature
+    estimate_a
+    analysis.estimate_a_classic
+    analysis.estimate_a_positive
 ```
 
 ## Other

notebooks/manual.ipynb

@@ -125,7 +125,7 @@
     "max_latitude = 48\n",
     "\n",
     "min_magnitude = 0.5\n",
-    "url = 'http://arclink.ethz.ch/fdsnws/event/1/query'\n",
+    "url = 'http://eida.ethz.ch/fdsnws/event/1/query'\n",
     "client = FDSNWSEventClient(url)\n",
     "df = client.get_events(\n",
     "    start_time=start_time,\n",

seismostats/__init__.py

@@ -2,9 +2,10 @@
 
 # flake8: noqa
 
-from seismostats.analysis.estimate_a import estimate_a
 # analysis
 from seismostats.analysis.estimate_beta import estimate_b, shi_bolt_confidence
+from seismostats.analysis.estimate_a import estimate_a
+
 # seismicity
 from seismostats.catalogs.catalog import Catalog, ForecastCatalog
 from seismostats.catalogs.rategrid import ForecastGRRateGrid, GRRateGrid

seismostats/analysis/__init__.py

@@ -3,7 +3,7 @@
 from seismostats.analysis.estimate_beta import (estimate_b_kijko_smit,
                                                 estimate_b_laplace,
                                                 estimate_b_positive,
-                                                estimate_b_tinti,
+                                                estimate_b_classic,
                                                 estimate_b_utsu,
                                                 estimate_b_weichert,
                                                 estimate_b_kijko_smit,

seismostats/analysis/estimate_a.py

@@ -1,91 +1,104 @@
 """This module contains functions for the estimation of a-value.
 """
 
-import numpy as np
 import warnings
+
+import numpy as np
+
 from seismostats.utils._config import get_option
 
 
 def estimate_a(magnitudes: np.ndarray,
-               times=None,
-               mc: float | None = None,
-               delta_m: float = None,
+               scaling_factor: float | None = None,
                m_ref: float | None = None,
+               mc: float | None = None,
                b_value: float | None = None,
-               scaling_factor: float | None = None,
+               delta_m: float = None,
+               times=None,
                method: str = "classic",
                ) -> float:
-    """Estimate the a-value of the Gutenberg-Richter (GR) law.
+    r"""Return the a-value of the Gutenberg-Richter (GR) law.
+
+    .. math::
+        N(m) = 10 ^ {a - b \cdot (m - m_{ref})}
+
+    where N(m) is the number of events with magnitude larger or equal to m
+    that occurred in the timeframe of the catalog.
 
     Args:
-        magnitudes: Magnitude sample
-        times:      Times of the events, in any format (datetime, float, etc.)
-        mc:         Completeness magnitude. If None, the lowest magnitude is
-                used as completeness magnitude.
-        delta_m:    Discretization of the magnitudes.
-        m_ref:      Reference magnitude for which the a-value is estimated. If
-                None, the a-value is estimated at mc.
-        b_value:    b-value of the Gutenberg-Richter distribution. Only needed
-                if m_ref is given.
-        scaling_factor:    Scaling factor. This should be chosen such that the
-                number of events observed can be normalized. For example:
-                Relative length of the time interval in which the events
-                occurred (relative to the time unit of interest, e.g., years).
+        magnitudes: vector of magnitudes, unsorted
+        scaling_factor:    scaling factor. If given, this is used to normalize
+                the number of observed events. For example:
+                Volume or area of the region considered
+                or length of the time interval, given in the unit of interest.
+        m_ref:      reference magnitude for which the a-value is estimated. If
+                None, m_ref is set to the smallest magnitude in the catalog.
+        mc:         completeness magnitude.
+                If given, magnitudes below mc are disregarded.
+                If None, the catalog is assumed to be complete and
+                mc is set to the smallest magnitude in the catalog.
+                This is only relevant when m_ref is not None.
+        b_value:    b-value of the Gutenberg-Richter law. Only relevant
+                when m_ref is not None.
+        delta_m:    discretization of magnitudes. default is no discretization
+        times:      vector of times of the events, in any format (datetime,
+                float, etc.) Only needed when method is "positive".
         method:     Method to estimate the a-value. Options are "classic" and
                 "positive".
 
+        Returns:
+        a: a-value of the Gutenberg-Richter law
+
     """
     if method == "classic":
         return estimate_a_classic(
-            magnitudes, mc, delta_m, m_ref, b_value, scaling_factor)
+            magnitudes, scaling_factor, m_ref, mc, b_value, delta_m)
     elif method == "positive":
         return estimate_a_positive(
             magnitudes,
             times,
-            mc,
             delta_m,
-            m_ref,
-            b_value,
-            scaling_factor,
+            scaling_factor=scaling_factor,
+            m_ref=m_ref,
+            mc=mc,
+            b_value=b_value,
             correction=False)
 
 
 def estimate_a_classic(magnitudes: np.ndarray,
-                       mc: float | None = None,
-                       delta_m: float = None,
+                       scaling_factor: float | None = None,
                        m_ref: float | None = None,
+                       mc: float | None = None,
                        b_value: float | None = None,
-                       scaling_factor: float | None = None,
+                       delta_m: float = None,
                        ) -> float:
-    """Estimate the a-value of the Gutenberg-Richter (GR) law.
+    r"""Return the a-value of the Gutenberg-Richter (GR) law.
 
-    N = 10 ** (a - b * (m_ref - mc)) (1)
+    .. math::
+        N(m) = 10 ^ {a - b \cdot (m - m_{ref})}
 
-    where N is the number of events with magnitude greater than m_ref, which
-    occurred in the timeframe of the catalogue.
-
-    If only the magnitudes are given, the a-value is estimated at the lowest
-    magnitude in the sample, with mc = min(magnitudes). Eq. (1) then simplifies
-    to N = 10**a.
+    where N(m) is the number of events with magnitude larger or equal to m
+    that occurred in the timeframe of the catalog.
 
     Args:
-        magnitudes: Magnitude sample
-        mc:         Completeness magnitude. If None, the lowest magnitude is
-                used as completeness magnitude.
-        delta_m:    Discretization of the magnitudes. This is needed solely to
-                avoid rounding errors. By default rounding errors are not
-                considered. This is adequate if the megnitudes are either
-                coninuous or do not contain rounding errors.
-        m_ref:      Reference magnitude for which the a-value is estimated. If
-                None, the a-value is estimated at mc.
-        b_value:    b-value of the Gutenberg-Richter distribution
-        scaling factor. This should be chosen such that the
-                number of events observed can be normalized. For example:
-                Relative length of the time interval in which the events
-                occurred (relative to the time unit of interest, e.g., years).
+        magnitudes: vector of magnitudes, unsorted
+        scaling_factor:    scaling factor. If given, this is used to normalize
+                the number of observed events. For example:
+                Volume or area of the region considered
+                or length of the time interval, given in the unit of interest.
+        m_ref:      reference magnitude for which the a-value is estimated. If
+                None, m_ref is set to the smallest magnitude in the catalog.
+        mc:         completeness magnitude.
+                If given, magnitudes below mc are disregarded.
+                If None, the catalog is assumed to be complete and
+                mc is set to the smallest magnitude in the catalog.
+                This is only relevant when m_ref is not None.
+        b_value:    b-value of the Gutenberg-Richter law. Only relevant
+                when m_ref is not None.
+        delta_m:    discretization of magnitudes. default is no discretization
 
     Returns:
-        a: a-value of the Gutenberg-Richter distribution
+        a: a-value of the Gutenberg-Richter law
     """
     if delta_m is None:
         delta_m = 0
@@ -118,35 +131,44 @@ def estimate_a_positive(
         magnitudes: np.ndarray,
         times: np.ndarray,
         delta_m: float,
-        mc: float | None = None,
         dmc: float | None = None,
+        scaling_factor: float | None = None,
         m_ref: float | None = None,
+        mc: float | None = None,
         b_value: float | None = None,
-        scaling_factor: float | None = None,
         correction: bool = False,
 ) -> float:
-    """Estimate the a-value of the Gutenberg-Richter (GR) law using only the
-    earthquakes whose magnitude m_i >= m_i-1 + dmc.
+    """Return the a-value of the Gutenberg-Richter (GR) law using only the
+    earthquakes with magnitude m_i >= m_i-1 + dmc.
+
+    Source:
+        Following the idea of positivity of
+        van der Elst 2021 (J Geophysical Research: Solid Earth, Vol 126, Issue
+        2).
+        Note: This is *not* a-positive as defined by van der Elst and Page 2023
+        (JGR: Solid Earth, Vol 128, Issue 10).
 
     Args:
-        magnitudes: Magnitude sample
-        times:      Times of the events. They should be sorted in ascending
-                order. Also, it is important that each time is scaled to the
-                time unit of interest (e.g., years). That is, in this case
-                each time should be a float and represent the time in years.
-        delta_m:    Discretization of the magnitudes.
-        dmc:        Minimum magnitude difference between consecutive events.
+        magnitudes: vector of magnitudes, unsorted
+        times:      vector of times of the events, in any format (datetime,
+                float, etc.)
+        delta_m:    discretization of the magnitudes.
+                default is no discretization
+        dmc:        minimum magnitude difference between consecutive events.
                 If None, the default value is delta_m.
-        mc:         Completeness magnitude. If None, the lowest magnitude is
-                used as completeness magnitude.
-        m_ref:      Reference magnitude for which the a-value is estimated. If
-                None, the a-value is estimated at mc.
-        b_value:    b-value of the Gutenberg-Richter distribution. Only needed
-                if m_ref is given.
-        scaling factor. This should be chosen such that the
-                number of events observed can be normalized. For example:
-                Relative length of the time interval in which the events
-                occurred (relative to the time unit of interest, e.g., years).
+        scaling_factor:    scaling factor. If given, this is used to normalize
+                the number of observed events. For example:
+                Volume or area of the region considered
+                or length of the time interval, given in the unit of interest.
+        m_ref:      reference magnitude for which the a-value is estimated. If
+                None, m_ref is set to the smallest magnitude in the catalog.
+        mc:         completeness magnitude.
+                If given, magnitudes below mc are disregarded.
+                If None, the catalog is assumed to be complete and
+                mc is set to the smallest magnitude in the catalog.
+                This is only relevant when m_ref is not None.
+        b_value:    b-value of the Gutenberg-Richter law. Only relevant
+                when m_ref is not None.
         correction: If True, the a-value is corrected for the bias introduced
                 by the observation period being larger than the time interval
                 between the first and last event. This is only relevant if the

seismostats/analysis/estimate_beta.py

@@ -41,7 +41,7 @@ def estimate_b(
     weights: list | None = None,
     b_parameter: str = "b_value",
     return_std: bool = False,
-    method="tinti",
+    method="classic",
     return_n: bool = False,
 ) -> float | tuple[float, float] | tuple[float, float, float]:
     """Return the maximum likelihood estimator for the Gutenberg-Richter
@@ -63,9 +63,9 @@ def estimate_b(
         method:     method to use for estimation of beta/b-value. Options
                 are:
 
-                - 'tinti',default, this is the is the classic estimator, see
-                  :func:`seismostats.analysis.estimate_b_tinti`
-                - 'positive' (this is b-positive, which applies the 'tinti'
+                - 'classic',default, this is the is the classic estimator, see
+                  :func:`seismostats.analysis.estimate_b_classic`
+                - 'positive' (this is b-positive, which applies the 'classic'
                   method to the positive differences, see
                   :func:`seismostats.analysis.estimate_b_positive`. To
                   achieve the effect of reduced STAI, the magnitudes must
@@ -98,8 +98,8 @@ def estimate_b(
                 "check if mc is chosen correctly"
             )
 
-    if method == "tinti":
-        return estimate_b_tinti(
+    if method == "classic":
+        return estimate_b_classic(
             magnitudes,
             mc=mc,
             delta_m=delta_m,
@@ -118,10 +118,10 @@ def estimate_b(
         )
 
     else:
-        raise ValueError("method must be either 'tinti' or 'positive'")
+        raise ValueError("method must be either 'classic' or 'positive'")
 
 
-def estimate_b_tinti(
+def estimate_b_classic(
     magnitudes: np.ndarray,
     mc: float,
     delta_m: float = 0,
@@ -301,7 +301,7 @@ def estimate_b_positive(
     # previous one. delta_m is added to avoid numerical errors
     mag_diffs = abs(mag_diffs[mag_diffs > dmc - delta_m / 2])
 
-    out = estimate_b_tinti(
+    out = estimate_b_classic(
         mag_diffs,
         mc=dmc,
         delta_m=delta_m,
@@ -375,7 +375,7 @@ def estimate_b_more_positive(
     # only take the values where the next earthquake is larger
     mag_diffs = abs(mag_diffs[mag_diffs > - delta_m / 2])
 
-    out = estimate_b_tinti(
+    out = estimate_b_classic(
         mag_diffs,
         mc=dmc,
         delta_m=delta_m,
@@ -484,7 +484,7 @@ def estimate_b_laplace(
     mag_diffs = abs(mag_diffs)
     mag_diffs = mag_diffs[mag_diffs > 0]
 
-    out = estimate_b_tinti(
+    out = estimate_b_classic(
         mag_diffs,
         mc=delta_m,
         delta_m=delta_m,

seismostats/analysis/estimate_mc.py

@@ -281,7 +281,7 @@ def mc_ks(
             warnings.warn("Both beta and b_method are given. Using beta.")
 
     if beta is None and b_method is None:
-        b_method = "tinti"
+        b_method = "classic"
 
     mcs_tested = []
     ks_ds = []