include dmc in b-more-pos

swiss-seismological-service · Jun 7, 2024 · ec1f1e4 · ec1f1e4
1 parent 18e45d4
commit ec1f1e4
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 6 deletions.
diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
@@ -318,6 +318,7 @@ def estimate_b_positive(
 def estimate_b_more_positive(
     magnitudes: np.ndarray,
     delta_m: float = 0,
+    dmc: float | None = None,
     b_parameter: str = "b_value",
     return_std: bool = False,
     return_n: bool = False,
@@ -352,12 +353,16 @@ def estimate_b_more_positive(
         std:    Shi and Bolt estimate of the beta/b-value estimate
         n:      number of events used for the estimation
     """
+
+    if dmc is None:
+        dmc = delta_m
+
     mag_diffs = np.zeros(len(magnitudes) - 1)
     for ii in range(len(magnitudes) - 1):
         for jj in range(ii + 1, len(magnitudes)):
             mag_diff_loop = magnitudes[jj] - magnitudes[ii]
             # print(mag_diff_loop, "diff loop")
-            if mag_diff_loop > 0:
+            if mag_diff_loop > dmc - delta_m / 2:
                 mag_diffs[ii] = mag_diff_loop
                 # print("take the value")
                 break
@@ -369,7 +374,7 @@ def estimate_b_more_positive(
 
     out = estimate_b_tinti(
         mag_diffs,
-        mc=delta_m,
+        mc=dmc,
         delta_m=delta_m,
         b_parameter=b_parameter,
         return_std=return_std,

diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py
@@ -199,10 +199,10 @@ def test_estimate_b_laplace(
 
 
 @pytest.mark.parametrize(
-    "n,b,mc,delta_m,b_parameter,precision",
+    "n,b,mc,delta_m,b_parameter,dmc,precision",
     [
-        (100000, 1.2 * np.log(10), 3, 0, "beta", 0.01),
-        (100000, 1, 3, 0.1, "b_value", 0.02),
+        (100000, 1.2 * np.log(10), 3, 0, "beta", None, 0.01),
+        (100000, 1, 3, 0.1, "b_value", 1, 0.02),
     ],
 )
 def test_estimate_b_more_positive(
@@ -211,13 +211,14 @@ def test_estimate_b_more_positive(
     mc: float,
     delta_m: float,
     b_parameter: str,
+    dmc: float,
     precision: float,
 ):
     mags = simulate_magnitudes_binned(
         n, b, mc, delta_m, b_parameter=b_parameter
     )
     b_estimate = estimate_b_more_positive(
-        mags, delta_m=delta_m, b_parameter=b_parameter
+        mags, delta_m=delta_m, dmc=dmc, b_parameter=b_parameter
     )
     assert abs(b - b_estimate) / b <= precision