shi and bolt uncertainty test function made

catalog_tools/__init__.py

@@ -6,7 +6,8 @@
     estimate_beta_elst,\
     estimate_beta_utsu,\
     estimate_beta_tinti,\
-    estimate_beta_laplace
+    estimate_beta_laplace,\
+    shi_bolt_confidence
 
 # download
 from catalog_tools.download.download_catalogs import\

catalog_tools/analysis/estimate_beta.py

@@ -121,7 +121,10 @@ def estimate_beta_laplace(
     return beta
 
 
-def shi_bolt_confidence(magnitudes: np.ndarray, b_value: float):
+def shi_bolt_confidence(
+        magnitudes: np.ndarray,
+        b_value: Optional[float] = None,
+        beta: Optional[float] = None):
     """ calculates the confidence limit according to shi and bolt 1982
 
     Source:
@@ -130,13 +133,22 @@ def shi_bolt_confidence(magnitudes: np.ndarray, b_value: float):
     Args:
         magnitudes: numpy array of magnitudes
         b_value:    b-value of the magnitudes
+        beta:       beta value (difference to b-value is factor of np.log(10)).
+                    -> provide either b_value or beta, not both
 
     Returns:
-        sig_b:  confidence limit of the b-value
+        sig_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
-    std_m = np.std(magnitudes, ddof=1) / np.sqrt(len(magnitudes))
-    sig_b = np.log(10) * b_value ** 2 * std_m
+    if b_value is not None:
+        std_m = np.std(magnitudes, ddof=1) / np.sqrt(len(magnitudes))
+        sig_b = np.log(10) * b_value ** 2 * std_m
+    elif beta is not None:
+        std_m = np.std(magnitudes, ddof=1) / np.sqrt(len(magnitudes))
+        sig_b = beta ** 2 * std_m
+    else:
+        print('input missing: b_value or beta')
+        sig_b = None
 
     return sig_b

catalog_tools/analysis/tests/test_estimate_beta.py

@@ -82,7 +82,41 @@ def test_estimate_beta_laplace(n: int, beta: float, mc: float, delta_m: float,
     beta_estimate = estimate_beta_laplace(mags, delta_m=delta_m)
     assert abs(beta - beta_estimate) / beta <= precision
 
-
-def test_shi_bolt_confidence():
-    shi_bolt_confidence(np.array([1, 1, 1, 1]), 1)
-    pass
+@pytest.mark.parametrize(
+    "magnitudes,b_value,std_b_value,std_beta",
+    [(np.array([0.20990507, 0.04077336, 0.27906596, 0.57406287, 0.64256544,
+       0.07293118, 0.58589873, 0.02678655, 0.27631233, 0.17682814]),
+      1, 0.16999880611649493, 0.39143671679062625),
+     (np.array([0.02637757, 0.06353823, 0.10257919, 0.54494906, 0.03928375,
+       0.08825028, 0.77713586, 0.54553981, 0.69315583, 0.06656642,
+       0.29035447, 0.2051877 , 0.30858087, 0.68896342, 0.03328782,
+       0.45016109, 0.40779409, 0.06788892, 0.02684032, 0.56140282,
+       0.29443359, 0.36328762, 0.17124489, 0.02154936, 0.36461541,
+       0.03613088, 0.15798366, 0.09111875, 0.16169287, 0.11986668,
+       0.10232035, 0.72695761, 0.19484174, 0.0459675 , 0.40514163,
+       0.08979514, 0.0442659 , 0.18672424, 0.21239088, 0.02287468,
+       0.1244267 , 0.04939361, 0.11232758, 0.02706083, 0.04275401,
+       0.02326945, 0.15048133, 0.50777581, 0.09583551, 0.40618488,
+       0.15595656, 0.09607254, 0.25576619, 0.01698973, 0.62755249,
+       0.31429311, 0.86575907, 0.37956298, 0.65648246, 0.0851286 ,
+       0.00850252, 0.22357953, 0.03295106, 0.08841752, 0.09657961,
+       0.54002676, 0.20335658, 0.23215333, 0.20120566, 0.60970099,
+       0.01128978, 0.31771308, 1.25246151, 0.02285632, 0.2687791 ,
+       0.1192099 , 0.06627574, 0.04301886, 0.24720467, 0.28518304,
+       0.04252851, 0.27818821, 0.08331663, 1.23090656, 0.1880176 ,
+       0.11314717, 0.01462853, 0.09256047, 0.4857446 , 0.08656431,
+       0.07022632, 0.32654491, 0.26047389, 0.20872121, 0.40157424,
+       0.02732529, 0.83884229, 0.4147758 , 0.07416183, 0.05636252]), 1.5,
+      0.13286469044352858, 0.3059322556005374)]
+)
+def test_shi_bolt_confidence(
+        magnitudes: np.ndarray,
+        b_value: float,
+        std_b_value: float,
+        std_beta: float):
+    precicion = 1e-10
+    beta = b_value * np.log(10)
+
+    assert shi_bolt_confidence(
+        magnitudes, b_value=b_value) - std_b_value < precicion
+    assert shi_bolt_confidence(magnitudes, beta=beta) - std_beta < precicion