Merge pull request #148 from swiss-seismological-service/documentatio…

…n/analysis

Documentation/analysis

docs/source/reference/analysis.md

@@ -0,0 +1,30 @@
+# Analysis
+
+```{eval-rst}
+.. currentmodule:: seismostats
+```
+
+## Estimating b-Values
+
+```{eval-rst}
+.. autosummary::
+    :toctree: api/
+
+    estimate_b
+    shi_bolt_confidence
+    analysis.estimate_b_tinti
+    analysis.estimate_b_positive
+    analysis.estimate_b_utsu
+    analysis.estimate_b_laplace
+    analysis.estimate_b_weichert
+    analysis.estimate_b_kijko_smit
+```
+
+## Estimating the Completeness Magnitude
+```{eval-rst}
+.. autosummary::
+    :toctree: api/
+
+    analysis.mc_ks
+    analysis.mc_max_curvature
+```

docs/source/reference/index.md

@@ -5,4 +5,5 @@
 
 catalog
 config
+analysis
 ```

seismostats/analysis/__init__.py

@@ -0,0 +1,9 @@
+# flake8: noqa
+
+from seismostats.analysis.estimate_beta import (estimate_b_kijko_smit,
+                                                estimate_b_laplace,
+                                                estimate_b_positive,
+                                                estimate_b_tinti,
+                                                estimate_b_utsu,
+                                                estimate_b_weichert)
+from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature

seismostats/analysis/estimate_beta.py

@@ -17,8 +17,10 @@ def estimate_b(
     method="tinti",
     return_n: bool = False,
 ) -> float | tuple[float, float] | tuple[float, float, float]:
-    """returns the maximum likelihood beta or b-value. Method depends on the
-    input parameter 'method'.
+    """Return the maximum likelihood estimator for the Gutenberg-Richter
+    b-value or beta, given an array of magnitudes
+    and a completeness magnitude.
+    Estimation method depends on the input parameter ``method``.
 
     Args:
         magnitudes: vector of magnitudes, unsorted, already cutoff (no
@@ -33,17 +35,24 @@ def estimate_b(
                 above) is returned
         method:     method to use for estimation of beta/b-value. Options
                 are:
-                    - 'tinti', 'utsu' (these are the classic estimators. 'tinti'
-                    is the recommended one, as it is more accurate. It is also
-                    the default method)
-                    - 'positive' (this is b-positive, which applies the 'tinti'
-                    method to the positive differences. To achieve the effect
-                    of reduced STAI, the magnitudes must be ordered in time)
-                    - 'laplace' (this is using the distribution of all
-                    differences, caution, can take long time to compute)
 
-        return_n:   if True the number of events used for the estimation is
-                returned. This is only relevant for the 'positive' method
+                - 'tinti', 'utsu' (these are the classic estimators, see
+                  :func:`seismostats.analysis.estimate_b_tinti`,
+                  :func:`seismostats.analysis.estimate_b_utsu`. 'tinti'
+                  is the recommended one, as it is more accurate. It is also
+                  the default method.)
+                - 'positive' (this is b-positive, which applies the 'tinti'
+                  method to the positive differences, see
+                  :func:`seismostats.analysis.estimate_b_positive`.
+                  To achieve the effect
+                  of reduced STAI, the magnitudes must be ordered in time.)
+                - 'laplace' (this is using the distribution of all
+                  differences, see
+                  :func:`seismostats.analysis.estimate_b_laplace`.
+                  Caution: it can take long time to compute.)
+
+        return_n:   If True, the number of events used for the estimation is
+                returned. This is only relevant for the 'positive' method.
 
     Returns:
         b:      maximum likelihood beta or b-value, depending on value of
@@ -116,11 +125,17 @@ def estimate_b_tinti(
     b_parameter: str = "b_value",
     return_std: bool = False,
 ) -> float | tuple[float, float]:
-    """returns the maximum likelihood beta
+    """Return the maximum likelihood b-value or beta for
+    an array of magnitudes and a completeness magnitude mc.
+    If the magnitudes are discretized, the discretization must be given in
+    ``delta_m``, so that the maximum likelihood estimator can be calculated
+    correctly.
+
+
     Source:
-        Aki 1965 (Bull. Earthquake research institute, vol 43, pp 237-239)
-        Tinti and Mulargia 1987 (Bulletin of the Seismological Society of
-            America, 77(6), 2125-2134.)
+        - Aki 1965 (Bull. Earthquake research institute, vol 43, pp 237-239)
+        - Tinti and Mulargia 1987 (Bulletin of the Seismological Society of
+          America, 77(6), 2125-2134.)
 
     Args:
         magnitudes: vector of magnitudes, unsorted, already cutoff (no
@@ -172,7 +187,8 @@ def estimate_b_utsu(
     b_parameter: str = "b_value",
     return_std: bool = False,
 ) -> float | tuple[float, float]:
-    """returns the maximum likelihood beta
+    """Return the maximum likelihood b-value or beta.
+
     Source:
         Utsu 1965 (Geophysical bulletin of the Hokkaido University, vol 13, pp
         99-103)
@@ -239,24 +255,26 @@ def estimate_b_positive(
     return_std: bool = False,
     return_n: bool = False,
 ) -> float | tuple[float, float] | tuple[float, float, float]:
-    """returns the b-value estimation using the positive differences of the
-    Magnitudes
+    """Return the b-value estimate calculated using the
+    positive differences between consecutive magnitudes.
 
     Source:
         Van der Elst 2021 (J Geophysical Research: Solid Earth, Vol 126, Issue
         2)
 
     Args:
-        magnitudes: vector of magnitudes differences, sorted in time (first
-                    entry is the earliest earthquake)
-        delta_m:    discretization of magnitudes. default is no discretization
+        magnitudes: array of magnitudes, sorted in time (first
+                    entry is the earliest earthquake).
+                    To achieve the effect
+                    of reduced STAI, the magnitudes must be ordered in time.
+        delta_m:    discretization of magnitudes. default is no discretization.
         b_parameter:either 'b-value', then the corresponding value  of the
                 Gutenberg-Richter law is returned, otherwise 'beta' from the
-                exponential distribution [p(M) = exp(-beta*(M-mc))]
+                exponential distribution [p(M) = exp(-beta*(M-mc))].
         return_std: if True the standard deviation of beta/b-value (see above)
-                is returned
+                is returned.
         return_n:   if True the number of events used for the estimation is
-                returned
+                returned.
 
     Returns:
         b:      maximum likelihood beta or b-value, depending on value of
@@ -293,9 +311,11 @@ def estimate_b_laplace(
     return_std: bool = False,
     return_n: bool = False,
 ) -> float | tuple[float, float]:
-    """returns the b-value estimation using the all the  differences of the
-    Magnitudes (this has a little less variance than the estimate_b_positive
-    method)
+    """Return the b-value estimate calculated using
+    all the  differences between magnitudes.
+    (This has a little less variance than the
+    :func:`seismostats.analysis.estimate_b_positive`
+    method.)
 
     Source:
         Van der Elst 2021 (J Geophysical Research: Solid Earth, Vol 126, Issue
@@ -339,10 +359,9 @@ def estimate_b_weichert(
     delta_m: float = 0.1,
     b_parameter: str = "b_value",
 ) -> tuple[float, float, float, float, float]:
-    """applies the Weichert (1980) algorithm for estimation of the
-    Gutenberg-Richter magnitude-frequency distribution parameters in
-    the case of unequal completeness periods for different magnitude
-    values.
+    """Return the b-value estimate calculated using the
+    Weichert (1980) algorithm, for the case of unequal
+    completeness periods for different magnitude values.
 
     Source:
         Weichert (1980), Estimation of the earthquake recurrence parameters
@@ -373,7 +392,7 @@ def estimate_b_weichert(
                     Gutenberg-Richter law is returned, otherwise 'beta'
                     from the exponential distribution [p(M) = exp(-beta*(M-mc))]
 
-    Returns:(
+    Returns:
         b_parameter: maximum likelihood point estimate of 'b-value' or 'beta'
         std_b_parameter: standard error of b_parameter
         rate_at_lmc: maximum likelihood point estimate of earthquake rate
@@ -557,23 +576,23 @@ def _weichert_objective_function(
     return np.abs(left - right)
 
 
-def estimate_b_kijko_smit_2012(
+def estimate_b_kijko_smit(
         magnitudes: np.ndarray,
         dates: list[np.datetime64],
         completeness_table: np.ndarray,
         last_year: int | float | None = None,
         delta_m: float = 0.1,
         b_parameter: str = 'b_value'
 ) -> tuple[float, float, float, float]:
-    """ applies the Kijko and Smit (2012) algorithm for estimation of the
-    Gutenberg-Richter magnitude-frequency distribution parameters in
-    the case of unequal completeness periods for different magnitude
-    values.
+    """Return the b-value estimate calculated using the
+    Kijko and Smit (2012) algorithm for the case of unequal
+    completeness periods for different magnitude values.
 
     Source:
-        Kijko A., Smit A. (2012), Extension of the Aki‐Utsu b‐Value Estimator
-        for Incomplete Catalogs, Bulletin of the Seismological Society of
-        America Vol 102, No. 3, pp. 1283–1287
+        Kijko, A. and Smit, A., 2012. Extension of the Aki‐Utsu b‐value
+        estimator for incomplete catalogs. Bulletin of the Seismological
+        Society of America, 102(3), pp.1283-1287.
+
 
     Args:
         magnitudes: vector of earthquake magnitudes
@@ -585,10 +604,10 @@ def estimate_b_kijko_smit_2012(
             the first column are considered detected. An example is given
             below:
 
-            np.array([[ 3.95, 1980],
-                      [ 4.95, 1920],
-                      [ 5.95, 1810],
-                      [ 6.95, 1520]])
+            >>> np.array([[ 3.95, 1980],
+            ...           [ 4.95, 1920],
+            ...           [ 5.95, 1810],
+            ...           [ 6.95, 1520]])
 
         last_year: last year of observation (the default is None, in which case
               it is set to the latest year in years).
@@ -597,7 +616,7 @@ def estimate_b_kijko_smit_2012(
                     Gutenberg-Richter law is returned, otherwise 'beta'
                     from the exponential distribution [p(M) = exp(-beta*(M-mc))]
 
-    Returns:(
+    Returns:
         b_parameter: maximum likelihood point estimate of 'b-value' or 'beta'
         std_b_parameter: standard error of b_parameter
         rate_at_lmc: maximum likelihood point estimate of earthquake rate
@@ -682,8 +701,8 @@ def shi_bolt_confidence(
     b: float | None = None,
     b_parameter: str = "b_value",
 ) -> float:
-    """calculates the confidence limit of the b_value or beta (depending on
-        which parameter is given) according to shi and bolt 1982
+    """Return the Shi and Bolt (1982) confidence limit of the b-value or
+    beta.
 
     Source:
         Shi and Bolt, BSSA, Vol. 72, No. 5, pp. 1677-1687, October 1982

seismostats/analysis/estimate_mc.py

@@ -2,15 +2,15 @@
 for the estimation of the completeness magnitude.
 """
 
+import warnings
+
 import numpy as np
 import pandas as pd
-import warnings
 
 from seismostats.analysis.estimate_beta import estimate_b
-from seismostats.utils.binning import get_fmd, bin_to_precision
-from seismostats.utils.simulate_distributions import (
-    simulated_magnitudes_binned,
-)
+from seismostats.utils.binning import bin_to_precision, get_fmd
+from seismostats.utils.simulate_distributions import \
+    simulated_magnitudes_binned
 
 
 def cdf_discrete_GR(
@@ -127,7 +127,7 @@ def ks_test_gr(
     )
 
     for ii in range(n):
-        simulated = simulated_all[n_sample * ii : n_sample * (ii + 1)]
+        simulated = simulated_all[n_sample * ii: n_sample * (ii + 1)]
         _, y_th = cdf_discrete_GR(simulated, mc=mc, delta_m=delta_m, beta=beta)
         _, y_emp = empirical_cdf(simulated)
 
@@ -155,9 +155,17 @@ def mc_ks(
     n: int = 10000,
 ) -> tuple[np.ndarray, list[float], np.ndarray, float | None, float | None]:
     """
-    Estimate the completeness magnitude (mc) for a given list of completeness
+    Return the completeness magnitude (mc) estimate
+    for a given list of completeness
     magnitudes using the K-S distance method.
 
+    Source:
+        - Clauset, A., Shalizi, C.R. and Newman, M.E., 2009. Power-law
+          distributions in empirical data. SIAM review, 51(4), pp.661-703.
+        - Mizrahi, L., Nandan, S. and Wiemer, S., 2021. The effect of
+          declustering on the size distribution of mainshocks. Seismological
+          Society of America, 92(4), pp.2333-2342.
+
     Args:
         sample:             Magnitudes to test
         delta_m:            Magnitude bins (sample has to be rounded to bins
@@ -278,8 +286,18 @@ def mc_max_curvature(
     correction_factor: float = 0.2,
 ) -> float:
     """
-    Estimate the completeness magnitude (mc) by maximum curvature (Wiemer and
-    Wyss 2000, Woessner and Wiemer 2005).
+    Return the completeness magnitude (mc) estimate
+    using the maximum curvature method.
+
+    Source:
+        - Wiemer, S. and Wyss, M., 2000. Minimum magnitude of completeness
+          in earthquake catalogs: Examples from Alaska, the western United
+          States, and Japan. Bulletin of the Seismological Society of America,
+          90(4), pp.859-869.
+        - Woessner, J. and Wiemer, S., 2005. Assessing the quality of earthquake
+          catalogues: Estimating the magnitude of completeness and its
+          uncertainty.
+          Bulletin of the Seismological Society of America, 95(2), pp.684-698.
     Args:
         sample:             Magnitudes to test
         delta_m:            Magnitude bins (sample has to be rounded to bins