Merge pull request #155 from swiss-seismological-service/feature/beta…

…2b_and_warnings

Feature/beta2b and warnings

seismostats/__init__.py

@@ -1,4 +1,5 @@
 """All functions that are to be used externally are initialized here"""
+
 # flake8: noqa
 
 # analysis

seismostats/analysis/estimate_beta.py

@@ -7,6 +7,32 @@
 import pandas as pd
 from scipy.optimize import minimize
 
+from seismostats.utils._config import get_option
+
+
+def beta_to_b_value(beta: float) -> float:
+    """converts the beta value to the b-value  of the Gutenberg-Richter law
+
+    Args:
+        beta: beta value
+
+    Returns:
+        b_value: corresponding b-value
+    """
+    return beta / np.log(10)
+
+
+def b_value_to_beta(b_value: float) -> float:
+    """converts the b-value to the beta value of the exponential distribution
+
+    Args:
+        b_value: b-value
+
+    Returns:
+        beta: corresponding beta value
+    """
+    return b_value * np.log(10)
+
 
 def estimate_b(
     magnitudes: np.ndarray,
@@ -37,23 +63,15 @@ def estimate_b(
         method:     method to use for estimation of beta/b-value. Options
                 are:
 
-                - '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.)
+                - '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'
                   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.
+                  :func:`seismostats.analysis.estimate_b_positive`. To
+                  achieve the effect of reduced STAI, the magnitudes must
+                  be ordered in time)
+        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
@@ -73,11 +91,12 @@ def estimate_b(
         np.min(magnitudes) >= mc
     ), "magnitudes below mc are present in the data"
     # test if lowest magnitude is much larger than mc
-    if np.min(magnitudes) - mc > delta_m / 2:
-        warnings.warn(
-            "no magnitudes in the lowest magnitude bin are present."
-            "check if mc is chosen correctly"
-        )
+    if get_option("warnings") is True:
+        if np.min(magnitudes) - mc > delta_m / 2:
+            warnings.warn(
+                "no magnitudes in the lowest magnitude bin are present."
+                "check if mc is chosen correctly"
+            )
 
     if method == "tinti":
         return estimate_b_tinti(
@@ -88,14 +107,7 @@ def estimate_b(
             b_parameter=b_parameter,
             return_std=return_std,
         )
-    elif method == "utsu":
-        return estimate_b_utsu(
-            magnitudes,
-            mc=mc,
-            delta_m=delta_m,
-            b_parameter=b_parameter,
-            return_std=return_std,
-        )
+
     elif method == "positive":
         return estimate_b_positive(
             magnitudes,
@@ -104,18 +116,9 @@ def estimate_b(
             return_std=return_std,
             return_n=return_n,
         )
-    elif method == "laplace":
-        return estimate_b_laplace(
-            magnitudes,
-            delta_m=delta_m,
-            b_parameter=b_parameter,
-            return_std=return_std,
-        )
 
     else:
-        raise ValueError(
-            "method must be either 'tinti', 'utsu', 'positive' or 'laplace'"
-        )
+        raise ValueError("method must be either 'tinti' or 'positive'")
 
 
 def estimate_b_tinti(
@@ -164,15 +167,14 @@ def estimate_b_tinti(
         beta = 1 / np.average(magnitudes - mc, weights=weights)
 
     if b_parameter == "b_value":
-        factor = 1 / np.log(10)
+        b = beta_to_b_value(beta)
     elif b_parameter == "beta":
-        factor = 1
+        b = beta
     else:
         raise ValueError(
             "please choose either 'b_value' or 'beta' as b_parameter"
         )
 
-    b = beta * factor
     std = shi_bolt_confidence(magnitudes, b=b, b_parameter=b_parameter)
 
     if return_std is True:
@@ -214,15 +216,14 @@ def estimate_b_utsu(
     beta = 1 / np.mean(magnitudes - mc + delta_m / 2)
 
     if b_parameter == "b_value":
-        factor = 1 / np.log(10)
+        b = beta_to_b_value(beta)
     elif b_parameter == "beta":
-        factor = 1
+        b = beta
     else:
         raise ValueError(
             "please choose either 'b_value' or 'beta' as b_parameter"
         )
 
-    b = beta * factor
     std = shi_bolt_confidence(magnitudes, b=b, b_parameter=b_parameter)
 
     if return_std is True:
@@ -282,6 +283,7 @@ def estimate_b_positive(
                 input variable 'b_parameter'. Note that the difference is just a
                 factor [b_value = beta * log10(e)]
         std:    Shi and Bolt estimate of the beta/b-value estimate
+        n:      number of events used for the estimation
     """
 
     mag_diffs = np.diff(magnitudes)
@@ -380,10 +382,10 @@ def estimate_b_weichert(
             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]])
 
         mag_max: maximum possible magnitude
         last_year: last year of observation (the default is None, in which case
@@ -415,16 +417,17 @@ def estimate_b_weichert(
         in np.arange(completeness_table[0, 0], mag_max + 0.001, delta_m)
         for i in np.unique(magnitudes)
     ], "magnitude bins not aligned with completeness edges"
-    if not np.all(magnitudes >= completeness_table[:, 0].min()):
-        warnings.warn(
-            "magnitudes below %.2f are not covered by the "
-            "completeness table and are discarded" % completeness_table[0, 0]
-        )
+    if get_option("warnings") is True:
+        if not np.all(magnitudes >= completeness_table[:, 0].min()):
+            warnings.warn(
+                "magnitudes below %.2f are not covered by the "
+                "completeness table and are discarded"
+                % completeness_table[0, 0]
+            )
     assert delta_m > 0, "delta_m cannot be zero"
     assert (
         b_parameter == "b_value" or b_parameter == "beta"
     ), "please choose either 'b_value' or 'beta' as b_parameter"
-    factor = 1 / np.log(10) if b_parameter == "b_value" else 1
 
     # convert datetime to integer calendar year
     years = np.array(dates).astype("datetime64[Y]").astype(int) + 1970
@@ -502,7 +505,6 @@ def estimate_b_weichert(
         tol=1e5 * np.finfo(float).eps,
     )
     beta = solution.x[0]
-    b_parameter = beta * factor
 
     # compute rate at lower magnitude of completeness bin
     weichert_multiplier = np.sum(
@@ -545,12 +547,15 @@ def estimate_b_weichert(
         * nominator
         / (denominator_term1 - denominator_term2)
     )
-    std_b_parameter = np.sqrt(var_beta) * factor
+    std_b = np.sqrt(var_beta)
+    if b_parameter == "b_value":
+        b = beta_to_b_value(beta)
+        std_b = beta_to_b_value(std_b)
 
     # compute uncertainty in rate at lower magnitude of completeness
     std_rate_at_lmc = rate_at_lmc / np.sqrt(complete_events.num.sum())
 
-    return b_parameter, std_b_parameter, rate_at_lmc, std_rate_at_lmc, a_val
+    return b, std_b, rate_at_lmc, std_rate_at_lmc, a_val
 
 
 def _weichert_objective_function(
@@ -589,12 +594,12 @@ def estimate_b_kijko_smit(
     delta_m: float = 0.1,
     b_parameter: str = "b_value",
 ) -> tuple[float, float, float, float]:
-    """Return the b-value estimate calculated using the
-    Kijko and Smit (2012) algorithm for 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. and Smit, A., 2012. Extension of the Aki‐Utsu b‐value
+        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.
 
@@ -647,16 +652,17 @@ def estimate_b_kijko_smit(
         )
         for i in np.unique(magnitudes)
     ], "magnitude bins not aligned with completeness edges"
-    if not np.all(magnitudes >= completeness_table[:, 0].min()):
-        warnings.warn(
-            "magnitudes below %.2f are not covered by the "
-            "completeness table and are discarded" % completeness_table[0, 0]
-        )
+    if get_option("warnings") is True:
+        if not np.all(magnitudes >= completeness_table[:, 0].min()):
+            warnings.warn(
+                "magnitudes below %.2f are not covered by the "
+                "completeness table and are discarded"
+                % completeness_table[0, 0]
+            )
     assert delta_m > 0, "delta_m cannot be zero"
     assert (
         b_parameter == "b_value" or b_parameter == "beta"
     ), "please choose either 'b_value' or 'beta' as b_parameter"
-    factor = 1 / np.log(10) if b_parameter == "b_value" else 1
 
     # convert datetime to integer calendar year
     years = np.array(dates).astype("datetime64[Y]").astype(int) + 1970
@@ -692,10 +698,13 @@ def estimate_b_kijko_smit(
         )
         estimator_terms.append((len(subcat_magnitudes) / ncomplete) / sub_beta)
     beta = 1 / np.sum(estimator_terms)
-    b_parameter = beta * factor
 
     # standard deviation of b/beta (Equation 8)
-    std_b_parameter = factor * beta / np.sqrt(ncomplete)
+    std_b = beta / np.sqrt(ncomplete)
+
+    if b_parameter == "b_value":
+        b = beta_to_b_value(beta)
+        std_b = beta_to_b_value(std_b)
 
     # get rate assuming Poisson process (Equation 10)
     denominator_rate = 0
@@ -710,7 +719,7 @@ def estimate_b_kijko_smit(
         completeness_magnitudes[0] + delta_m * 0.5
     )
 
-    return b_parameter, std_b_parameter, rate_at_lmc, a_val
+    return b, std_b, rate_at_lmc, a_val
 
 
 def shi_bolt_confidence(
@@ -730,20 +739,18 @@ def shi_bolt_confidence(
         b_parameter:either either 'b_value' or 'beta'
 
     Returns:
-        std:  confidence limit of the b-value/beta value (depending on input)
+        std_b:  confidence limit of the b-value/beta value (depending on input)
     """
     # standard deviation in Shi and Bolt is calculated with 1/(N*(N-1)), which
     # is by a factor of sqrt(N) different to the std(x, ddof=1) estimator
-    if b_parameter == "b_value":
-        factor = 1 / np.log(10)
-    elif b_parameter == "beta":
-        factor = 1
-    else:
-        raise ValueError(
-            "please choose either 'b_value' or 'beta' as b_parameter"
-        )
+    assert (
+        b_parameter == "b_value" or b_parameter == "beta"
+    ), "please choose either 'b_value' or 'beta' as b_parameter"
 
-    std = np.std(magnitudes, ddof=1) / np.sqrt(len(magnitudes))
-    std = 1 / factor * b**2 * std
+    std_b = (
+        np.log(10) * b**2 * np.std(magnitudes) / np.sqrt(len(magnitudes) - 1)
+    )
+    if b_parameter == "beta":
+        std_b = (std_b) / np.log(10)
 
-    return std
+    return std_b