Feature/b val uncertainty

.gitignore

@@ -133,3 +133,6 @@ dmypy.json
 
 # Swap files
 .swp
+
+# Finder files
+.DS_Store

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

@@ -4,8 +4,13 @@
 from typing import Optional
 
 
-def estimate_beta_tinti(magnitudes: np.ndarray, mc: float, delta_m: float = 0,
-                        weights: Optional[list] = None) -> float:
+def estimate_beta_tinti(magnitudes: np.ndarray,
+                        mc: float,
+                        delta_m: float = 0,
+                        weights: Optional[list] = None,
+                        gutenberg: bool = False,
+                        error: bool = False
+                        ) -> (float, float):
     """ returns the maximum likelihood beta
     Source:
         Aki 1965 (Bull. Earthquake research institute, vol 43, pp 237-239)
@@ -18,9 +23,16 @@ def estimate_beta_tinti(magnitudes: np.ndarray, mc: float, delta_m: float = 0,
         mc:         completeness magnitude
         delta_m:    discretization of magnitudes. default is no discretization
         weights:    weights of each magnitude can be specified here
+        gutenberg:  if True the b-value of the Gutenberg-Richter law is
+                    returned, otherwise the beta value of the exponential
+                    distribution [p(M) = exp(-beta*(M-mc))] is returned
+        error:      if True the error of beta/b-value (see above) is returned
 
     Returns:
-        beta:       maximum likelihood beta (b_value = beta * log10(e))
+        b:          maximum likelihood beta or b-value, depending on value of
+                    input variable 'gutenberg'. Note that the difference
+                    is just a factor [b_value = beta * log10(e)]
+        std_b:      Shi and Bolt estimate of the beta/b-value estimate
     """
 
     if delta_m > 0:
@@ -29,7 +41,18 @@ def estimate_beta_tinti(magnitudes: np.ndarray, mc: float, delta_m: float = 0,
     else:
         beta = 1 / np.average(magnitudes - mc, weights=weights)
 
-    return beta
+    if gutenberg is True:
+        factor = 1 / np.log(10)
+    else:
+        factor = 1
+
+    if error is True:
+        std_b = shi_bolt_confidence(magnitudes, beta=beta) * factor
+        b = beta * factor
+        return b, std_b
+    else:
+        b = beta * factor
+        return b
 
 
 def estimate_beta_utsu(magnitudes: np.ndarray, mc: float, delta_m: float = 0
@@ -119,3 +142,38 @@ def estimate_beta_laplace(
     mag_diffs = mag_diffs[mag_diffs > 0]
     beta = estimate_beta_tinti(mag_diffs, mc=delta_m, delta_m=delta_m)
     return beta
+
+
+def shi_bolt_confidence(
+        magnitudes: np.ndarray,
+        b_value: Optional[float] = None,
+        beta: Optional[float] = None
+) -> float:
+    """ calculates the confidence limit of the b_value or beta (depending on
+        which parameter is given) according to shi and bolt 1982
+
+    Source:
+        Shi and Bolt, BSSA, Vol. 72, No. 5, pp. 1677-1687, October 1982
+
+    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/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_value is not None:
+        std_m = np.std(magnitudes, ddof=1) / np.sqrt(len(magnitudes))
+        std_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))
+        std_b = beta ** 2 * std_m
+    else:
+        print('input missing: b_value or beta')
+        std_b = None
+
+    return std_b

catalog_tools/analysis/tests/test_estimate_beta.py

@@ -11,7 +11,8 @@
     estimate_beta_utsu,\
     estimate_beta_elst,\
     estimate_beta_laplace,\
-    differences
+    differences,\
+    shi_bolt_confidence
 
 
 def simulate_magnitudes_w_offset(n: int, beta: float, mc: float,
@@ -80,3 +81,43 @@ def test_estimate_beta_laplace(n: int, beta: float, mc: float, delta_m: float,
     mags = simulate_magnitudes_w_offset(n, beta, mc, delta_m)
     beta_estimate = estimate_beta_laplace(mags, delta_m=delta_m)
     assert abs(beta - beta_estimate) / beta <= precision
+
+
+@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):
+    precision = 1e-10
+    beta = b_value * np.log(10)
+
+    assert shi_bolt_confidence(
+        magnitudes, b_value=b_value) - std_b_value < precision
+    assert shi_bolt_confidence(magnitudes, beta=beta) - std_beta < precision

catalog_tools/download/download_catalogs.py

@@ -179,12 +179,6 @@ def download_catalog_1(
     data = response.read()
     df = pd.read_csv(BytesIO(data), delimiter=delimiter)
 
-    # for saving the byte file:
-    # with open(
-    #         '/Users/aron/polybox/Projects/catalog-tools/catalog_tools/
-    #         download/tests/data/catalog_scedc.bin',
-    #         'wb') as f:
-    #     f.write(data)
     return df
 
 
@@ -241,8 +235,6 @@ def prepare_scedc_catalog(
         delta_m: magnitude bin size to be applied.
         only_earthquakes: if True, only
             events of event_type earthquake are kept.
-        convert_to_mw: if True, local magnitudes are converted to Mw
-            using Edwards et al.
 
     Returns:
         the catalog as a DataFrame
@@ -291,8 +283,8 @@ def download_catalog(
         min_longitude:  minimum longitude of catalog.
         max_longitude:  maximum longitude of catalog.
         min_magnitude:  minimum magnitude of catalog.
-        delta_m:        magnitude bin size. if >0, then events of
-            magnitude >= (min_magnitude - delta_m/2) will be downloaded.
+
+    TODO: include delta_m in parameters?
 
     Returns:
         The catalog as a pandas DataFrame.
@@ -385,7 +377,6 @@ def download_catalog(
         evs.append(pd.Series([time, lat, lon, depth, mag, mag_type]))
 
     cat = pd.DataFrame(evs)
-    # cat.to_csv('/Users/aron/polybox/Projects/catalog-tools/catalog_tools/download/tests/data/catalog.csv')
 
     cat.columns = [
         'time', 'latitude', 'longitude', 'depth', 'magnitude', 'mag_type']