merge with main

seismostats/analysis/bvalue/base.py

@@ -7,8 +7,8 @@
 
 from seismostats.analysis.bvalue.utils import (b_value_to_beta,
                                                shi_bolt_confidence)
-from seismostats.utils.binning import binning_test
 from seismostats.utils._config import get_option
+from seismostats.utils.binning import binning_test
 
 
 class BValueEstimator(ABC):
@@ -72,7 +72,16 @@ def std(self):
         Shi and Bolt estimate of the beta/b-value estimate.
         '''
         assert self.__b_value is not None, 'Please run the estimator first.'
+        if get_option('warnings') is True:
+            if self.weights is not None:
+                warnings.warn(
+                    'Shi and Bolt confidence with weights considers the '
+                    'magnitudes as '
+                    'having length {}, the sum of relevant weights.'.format(
+                        np.sum(self.weights))
+                )
         return shi_bolt_confidence(self.magnitudes,
+                                   weights=self.weights,
                                    b=self.__b_value,
                                    b_parameter=self.__b_parameter)
 
@@ -102,6 +111,12 @@ def _sanity_checks(self):
         )
         "Magnitudes are not binned correctly."
 
+        if self.weights is not None:
+            assert len(self.magnitudes) == len(self.weights), (
+                "The number of magnitudes and weights must be equal."
+            )
+            assert np.all(self.weights >= 0), "Weights must be nonnegative."
+
         # test if lowest magnitude is much larger than mc
         if get_option("warnings") is True:
             if np.min(self.magnitudes) - self.mc > self.delta_m / 2:

seismostats/analysis/bvalue/more_positive.py

@@ -23,7 +23,7 @@ class BMorePositiveBValueEstimator(BValueEstimator):
                     value are not considered). If None, it is set to delta_m.
     """
 
-    weights_supported = False
+    weights_supported = True
 
     def __init__(self, dmc: float | None = None, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -38,17 +38,25 @@ def __init__(self, dmc: float | None = None, *args, **kwargs):
     def _estimate(self):
 
         mag_diffs = - np.ones(len(self.magnitudes) - 1) * self.delta_m
+        if self.weights is not None:
+            # weight vector of same length as mag diffs
+            weights = - np.ones(len(self.magnitudes) - 1) * self.delta_m
         for ii in range(len(self.magnitudes) - 1):
             for jj in range(ii + 1, len(self.magnitudes)):
                 mag_diff_loop = self.magnitudes[jj] - self.magnitudes[ii]
                 if mag_diff_loop > self.dmc - self.delta_m / 2:
                     mag_diffs[ii] = mag_diff_loop
+                    if self.weights is not None:
+                        # use weight of second earthquake of a difference
+                        weights[ii] = self.weights[jj]
                     break
 
         # only take the values where the next earthquake is larger
+        if self.weights is not None:
+            self.weights = weights[mag_diffs > - self.delta_m / 2]
         self.magnitudes = abs(mag_diffs[mag_diffs > - self.delta_m / 2])
 
         classic_estimator = ClassicBValueEstimator(mc=self.dmc,
                                                    delta_m=self.delta_m)
 
-        return classic_estimator(self.magnitudes)
+        return classic_estimator(self.magnitudes, weights=self.weights)

seismostats/analysis/bvalue/positive.py

@@ -21,7 +21,7 @@ class BPositiveBValueEstimator(BValueEstimator):
                 value are not considered). If None, the cutoff is set to delta_m
     '''
 
-    weights_supported = False
+    weights_supported = True
 
     def __init__(self, dmc: float | None = None, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -37,10 +37,15 @@ def _estimate(self):
         self.magnitudes = np.diff(self.magnitudes)
         # only take the values where the next earthquake is d_mc larger than the
         # previous one. delta_m is added to avoid numerical errors
+        if self.weights is not None:
+            # use weight of second earthquake of a difference
+            self.weights = self.weights[1:]
+            self.weights = self.weights[self.magnitudes
+                                        > self.dmc - self.delta_m / 2]
         self.magnitudes = abs(
             self.magnitudes[self.magnitudes > self.dmc - self.delta_m / 2])
 
         classic_estimator = ClassicBValueEstimator(mc=self.dmc,
                                                    delta_m=self.delta_m)
 
-        return classic_estimator(self.magnitudes)
+        return classic_estimator(self.magnitudes, self.weights)

seismostats/analysis/bvalue/tests/test_bvalues.py

@@ -42,10 +42,20 @@ def test_estimate_b_classic(
     delta_m: float,
 ):
     mags = bin_to_precision(mags, delta_m)
+    mags_extended = np.concatenate([mags, mags + delta_m])
+    weights = np.ones(len(mags))
+    weights_extended = np.concatenate([weights, 0 * weights])
+
     estimator = ClassicBValueEstimator(mc=mc, delta_m=delta_m)
     b_estimate = estimator(mags)
+    b_estimate_weighted = estimator(mags, weights=weights)
+    b_estimate_half_weighted = estimator(mags, weights=weights * 0.5)
+    b_estimate_extended = estimator(mags_extended, weights=weights_extended)
 
     assert_almost_equal(b_estimate, b_est_correct)
+    assert_almost_equal(b_estimate, b_estimate_weighted)
+    assert_almost_equal(b_estimate, b_estimate_half_weighted)
+    assert_almost_equal(b_estimate, b_estimate_extended)
 
 
 @pytest.mark.parametrize(
@@ -63,9 +73,20 @@ def test_estimate_b_utsu(
     delta_m: float,
 ):
     mags = bin_to_precision(mags, delta_m)
+    mags_extended = np.concatenate([mags, mags + delta_m])
+    weights = np.ones(len(mags))
+    weights_extended = np.concatenate([weights, 0 * weights])
+
     estimator = UtsuBValueEstimator(mc=mc, delta_m=delta_m)
     b_estimate = estimator(mags)
+    b_estimate_weighted = estimator(mags, weights=weights)
+    b_estimate_half_weighted = estimator(mags, weights=weights * 0.5)
+    b_estimate_extended = estimator(mags_extended, weights=weights_extended)
+
     assert_almost_equal(b_estimate, b_est_correct)
+    assert_almost_equal(b_estimate, b_estimate_weighted)
+    assert_almost_equal(b_estimate, b_estimate_half_weighted)
+    assert_almost_equal(b_estimate, b_estimate_extended)
 
 
 @pytest.mark.parametrize(
@@ -84,9 +105,21 @@ def test_estimate_b_positive(
     dmc: float,
 ):
     mags = bin_to_precision(mags, delta_m)
+    mags = mags[mags >= mc - delta_m / 2]
+    mags_extended = np.concatenate([mags, mags + delta_m])
+    weights = np.ones(len(mags))
+    weights_extended = np.concatenate([weights, 0 * weights])
+
     estimator = BPositiveBValueEstimator(mc=mc, delta_m=delta_m, dmc=dmc)
     b_estimate = estimator(mags)
+    b_estimate_weighted = estimator(mags, weights=weights)
+    b_estimate_half_weighted = estimator(mags, weights=weights * 0.5)
+    b_estimate_extended = estimator(mags_extended, weights=weights_extended)
+
     assert_almost_equal(b_estimate, b_est_correct)
+    assert_almost_equal(b_estimate, b_estimate_weighted)
+    assert_almost_equal(b_estimate, b_estimate_half_weighted)
+    assert_almost_equal(b_estimate, b_estimate_extended)
 
 
 @pytest.mark.parametrize(
@@ -105,6 +138,14 @@ def test_estimate_b_more_positive(
     dmc: float
 ):
     mags = bin_to_precision(mags, delta_m)
+    mags = mags[mags >= mc - delta_m / 2]
+    weights = np.ones(len(mags))
+
     estimator = BMorePositiveBValueEstimator(mc=mc, delta_m=delta_m, dmc=dmc)
     b_estimate = estimator(mags)
+    b_estimate_weighted = estimator(mags, weights=weights)
+    b_estimate_half_weighted = estimator(mags, weights=weights * 0.5)
+
     assert_almost_equal(b_estimate, b_est_correct)
+    assert_almost_equal(b_estimate, b_estimate_weighted)
+    assert_almost_equal(b_estimate, b_estimate_half_weighted)

seismostats/analysis/bvalue/tests/test_utils.py

@@ -1,8 +1,8 @@
-from numpy.testing import assert_almost_equal
 from datetime import datetime
 
 import numpy as np
 import pytest
+from numpy.testing import assert_almost_equal
 
 from seismostats.analysis.bvalue.tests.test_bvalues import magnitudes
 from seismostats.analysis.bvalue.utils import (b_value_to_beta,
@@ -55,5 +55,14 @@ def test_make_more_incomplete():
 )
 def test_shi_bolt_confidence(
         std: float, mags: np.ndarray, b: float, b_parameter: str):
-    assert_almost_equal(
-        shi_bolt_confidence(mags, b=b, b_parameter=b_parameter), std)
+    weights = np.ones(len(mags))
+
+    conf = shi_bolt_confidence(mags, b=b, b_parameter=b_parameter)
+    conf_weighted = shi_bolt_confidence(mags, weights=weights, b=b,
+                                        b_parameter=b_parameter)
+    conf_half_weighted = shi_bolt_confidence(
+        mags, weights=weights * 0.5, b=b, b_parameter=b_parameter)
+
+    assert_almost_equal(conf_weighted, std)
+    assert_almost_equal(conf, std)
+    assert (conf_half_weighted > conf)

seismostats/analysis/bvalue/utils.py

@@ -27,6 +27,7 @@ def b_value_to_beta(b_value: float) -> float:
 
 def shi_bolt_confidence(
     magnitudes: np.ndarray,
+    weights: np.ndarray | None = None,
     b: float | None = None,
     b_parameter: str = "b_value",
 ) -> float:
@@ -38,6 +39,7 @@ def shi_bolt_confidence(
 
     Args:
         magnitudes: numpy array of magnitudes
+        weights:    numpy array of weights for each magnitude
         b:          known or estimated b-value/beta of the magnitudes
         b_parameter:either either 'b_value' or 'beta'
 
@@ -50,8 +52,15 @@ def shi_bolt_confidence(
         b_parameter == "b_value" or b_parameter == "beta"
     ), "please choose either 'b_value' or 'beta' as b_parameter"
 
+    std_mags = np.sqrt(np.average(np.square(
+        magnitudes - np.average(magnitudes, weights=weights)), weights=weights))
+    if weights is None:
+        len_mags = len(magnitudes)
+    else:
+        len_mags = np.sum(weights)
+
     std_b = (
-        np.log(10) * b**2 * np.std(magnitudes) / np.sqrt(len(magnitudes) - 1)
+        np.log(10) * b**2 * std_mags / np.sqrt(len_mags - 1)
     )
     if b_parameter == "beta":
         std_b = (std_b) / np.log(10)

seismostats/analysis/bvalue/utsu.py

@@ -6,7 +6,7 @@
 
 class UtsuBValueEstimator(BValueEstimator):
 
-    weights_supported = False
+    weights_supported = True
 
     def __init__(self, *args, **kwargs):
         """Return the maximum likelihood b-value or beta.
@@ -24,5 +24,6 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
     def _estimate(self):
-        beta = 1 / np.mean(self.magnitudes - self.mc + self.delta_m / 2)
+        beta = 1 / np.average(self.magnitudes - self.mc
+                              + self.delta_m / 2, weights=self.weights)
         return beta_to_b_value(beta)