added plotting functionality

seismostats/analysis/b_significant.py

@@ -136,24 +136,27 @@ def transform_n(
     return b_transformed
 
 
-def b_series(
+def bs_from_partitioning(
     list_magnitudes: list[np.ndarray],
     list_times: list[np.ndarray[np.datetime64]],
     delta_m: float,
     mc: float | np.ndarray,
     b_method: BValueEstimator = ClassicBValueEstimator,
-
-
+    **kwargs,
 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
-    """ Estimate the series of b-values
+    """ Estimate the series of b-values from a list of subsets of magnitudes and
+    times.
 
     Args:
-        list_mags:  list of arrays of magnitudes
-        list_times: list of arrays of times
+        list_mags:  list of arrays of magnitudes. Example of how to provide the
+            data: [np.array([1, 2, 3], np.array([2, 3, 4])]. From each array
+            within the list, a b-value is estimated.
+        list_times: list of arrays of times, in the same order as the magnitudes
         delta_m:    discretization of magnitudes
         mc:         completeness magnitude (can be a vector of same
                 length as list_magnitudes)
         b_method:   method to estimate the b-value
+        **kwargs:   additional arguments to the b-value estimation method
 
     Returns:
         b_values:   series of b-values, each one is estimated from the
@@ -163,7 +166,7 @@ def b_series(
     """
 
     b_values = np.zeros(len(list_magnitudes))
-    std_b = np.zeros(len(list_magnitudes))
+    std_bs = np.zeros(len(list_magnitudes))
     n_ms = np.zeros(len(list_magnitudes))
     if isinstance(mc, (float, int)):
         mc = np.ones(len(list_magnitudes)) * mc
@@ -177,12 +180,12 @@ def b_series(
         mags_loop = mags_loop[idx_sorted]
         times_loop = times_loop[idx_sorted]
 
-        estimator.calculate(mags_loop, mc=mc[ii], delta_m=delta_m)
+        estimator.calculate(mags_loop, mc=mc[ii], delta_m=delta_m, **kwargs)
         b_values[ii] = estimator.b_value
-        std_b[ii] = estimator.std
+        std_bs[ii] = estimator.std
         n_ms[ii] = estimator.n
 
-    return b_values, std_b, n_ms.astype(int)
+    return b_values, std_bs, n_ms.astype(int)
 
 
 def cut_constant_idx(
@@ -228,6 +231,7 @@ def mac_1D_constant_nm(
         n_m: int,
         min_num: int = 10,
         b_method: BValueEstimator = ClassicBValueEstimator,
+        **kwargs,
 ) -> tuple[float, float, float, np.ndarray, np.ndarray]:
     """
     This function estimates the mean autocorrelation for the one-dimensional
@@ -241,6 +245,10 @@ def mac_1D_constant_nm(
     As a lower limit, no less than 25 subsamples (which can be estimated by
     len(mags) / n_m) should be used.
 
+    In order to plot the corresponding series of b-values, use the function
+    plot_b_constant_mn (from seismostats.plots.statistical  import
+    plot_b_constant_mn) with the same parameters as used here.
+
     Args:
         mags:   magnitudes of the events. They are assumed to be order along the
             dimension of interest (e.g. time or depth)
@@ -252,11 +260,10 @@ def mac_1D_constant_nm(
         n_m:   number of magnitudes in each partition
         min_num:    minimum number of events in a partition
         b_method:   method to estimate the b-values
-        return_nm:  if True, the mean number of events per b-value estimate is
-            returned
+        **kwargs:   additional arguments to the b-value estimation method
 
     Returns:
-        mac:        mean autocorrelation.
+        mac:        mean autocorrelation
         mu_mac:     expected mean autocorrelation und H0
         std_mac:    standard deviation of the mean autocorrelation under H0
                 (i.e. constant b-value). Here, the conservatice estimate
@@ -293,52 +300,49 @@ def mac_1D_constant_nm(
     if len(mags) != len(times):
         raise ValueError("mags and times must have the same length")
 
-    # estimate a and b values for n realizations
+    # estimate a and b values for n_m realizations
     ac_1D = np.zeros(n_m)
     n = np.zeros(n_m)
     n_p = np.zeros(n_m)
     n_ms = np.zeros(n_m)
 
     for ii in range(n_m):
         # partition data
-        idx_left, tile_magnitudes = cut_constant_idx(
+        idx, list_magnitudes = cut_constant_idx(
             mags, n_m, offset=ii
         )
-        tile_times = np.array_split(times, idx_left)
-        tile_mc = np.array_split(mc, idx_left)
-        for jj, mc_loop in enumerate(tile_mc):
-            tile_mc[jj] = float(max(mc_loop))
+        list_times = np.array_split(times, idx)
+        list_mc = np.array_split(mc, idx)
+        for jj, mc_loop in enumerate(list_mc):
+            list_mc[jj] = float(max(mc_loop))
 
         # make sure that data at the edges is not included if not enough
         # samples
-        if len(tile_magnitudes[-1]) < n_m:
-            tile_magnitudes.pop(-1)
-            tile_times.pop(-1)
-            tile_mc.pop(-1)
-
-        if len(tile_magnitudes[0]) < n_m:
-            tile_magnitudes.pop(0)
-            tile_times.pop(0)
-            tile_mc.pop(0)
+        if len(list_magnitudes[-1]) < n_m:
+            list_magnitudes.pop(-1)
+            list_times.pop(-1)
+            list_mc.pop(-1)
+        if len(list_magnitudes[0]) < n_m:
+            list_magnitudes.pop(0)
+            list_times.pop(0)
+            list_mc.pop(0)
 
         # estimate b-values
-        b_vec, _, n_m_loop = b_series(
-            tile_magnitudes, tile_times, delta_m,
-            tile_mc, b_method=b_method)
+        b_vec, _, n_m_loop = bs_from_partitioning(
+            list_magnitudes, list_times, delta_m,
+            list_mc, b_method=b_method, **kwargs)
         b_vec[n_m_loop < min_num] = np.nan
 
         # estimate average events per b-value estimate
         n_ms[ii] = np.mean(n_m_loop[n_m_loop >= min_num])
-
-        # estimate autocorrelation (1D, not considering nan)
+        # estimate autocorrelation (1D)
         ac_1D[ii], n[ii], n_p[ii], = est_morans_i(b_vec)
 
+    # estimate mean and (conservative) standard deviation of the
+    # autocorrelation under H0
     mac = np.nanmean(ac_1D)
     mean_n = np.nanmean(n)
     mean_np = np.nanmean(n_p)
-
-    # estimate mean and (conservative )standard deviation of the
-    # autocorrelation under H0
     mu_mac = -1 / mean_n
     std_mac = 1 / np.sqrt(mean_np)
 

seismostats/analysis/tests/test_b_significant.py

@@ -5,7 +5,7 @@
 
 from seismostats.analysis.b_significant import (
     est_morans_i,
-    b_series,
+    bs_from_partitioning,
     cut_constant_idx,
     transform_n,
 )
@@ -72,7 +72,7 @@ def test_transform_n():
         transform_n(b_est, b_true, n1, np.array([1, 2]))
 
 
-def test_b_series():
+def test_bs_from_partitioning():
     list_magnitudes = [np.array([1, 2, 3, 4, 5]),
                        np.array([1, 2, 3, 4, 5]),
                        np.array([1, 2, 3, 4, 5])]
@@ -87,7 +87,8 @@ def test_b_series():
                             datetime(2021, 1, 15)])]
     delta_m = 1
     mc = 1
-    b_values, std_b, n_ms = b_series(list_magnitudes, list_times, delta_m, mc)
+    b_values, std_b, n_ms = bs_from_partitioning(
+        list_magnitudes, list_times, delta_m, mc)
     assert_almost_equal(b_values, np.array(
         [0.17609125905568124, 0.17609125905568124, 0.17609125905568124]))
     assert_almost_equal(std_b, np.array([0.05048661905780697,

seismostats/plots/statistical.py

@@ -1,6 +1,8 @@
 # standard
 import matplotlib.pyplot as plt
 import numpy as np
+from typing import Literal
+
 # statistical
 from scipy.stats import norm
 
@@ -71,3 +73,137 @@ def plot_mc_vs_b(
         ax.legend()
 
     return ax
+
+
+def plot_b_series_constant_nm(
+        mags: np.ndarray,
+        delta_m: float,
+        mc: np.ndarray,
+        times: np.ndarray,
+        n_m: int,
+        min_num: float = 2,
+        b_method: BValueEstimator = ClassicBValueEstimator,
+        plot_technique: Literal['left', 'midpoint', 'right'] = 'right',
+        x_variable: np.ndarray | None = None,
+        confidence: float = 0.95,
+        ax: plt.Axes | None = None,
+        color: str = "blue",
+        label: str | None = None,
+        **kwargs,
+) -> plt.Axes:
+    """
+    Plots the b-values estimated from a running window of n_m magnitudes.
+
+    Args:
+        mags:   magnitudes of the events. If x_variable is None, the magnitudes
+            are assumed to be sorted in the dimension of interest.
+        delta_m:    magnitude bin width
+        mc:     completeness magnitude. If a single value is provided, it is
+            used for all magnitudes. Otherwise, the individual completeness of
+            each magnitude can be provided.
+        times:  times of the events
+        n_m:   number of magnitudes in each partition
+        min_num:    minimum number of events from which a b-value is estimated.
+            If the number of events is smaller, the b-value is set to np.nan
+        b_method:   method to estimate the b-values
+        plot_technique:    technique where to plot the b-values with respect to
+            the x-variable. Options are 'left', 'midpoint', 'right'. If set to
+            'right' (default), the b-value is plotted at the right edge. For
+            time series, this is the most common choice, as it avoids optical
+            illusions of the b-value predicting future seismicity.
+        x_variable: values of the dimension of interest, along which the
+            b-values should be plotted. It should be a 1D array with the same
+            length as the magnitudes, e.g., the time of the events. If None,
+            the b-values are plotted against the event index.
+        confidence:    confidence interval that should be plotted. Default
+            is 0.95 (i.e., the 95% confidence interval is plotted)
+        ax:    axis where the plot should be plotted
+        color: color of the data
+
+    Returns:
+        ax that was plotted on
+    """
+
+    if isinstance(mc, (float, int)):
+        if min(mags) < mc:
+            raise ValueError("The completeness magnitude is larger than the "
+                             "smallest magnitude")
+        mc = np.ones(len(mags)) * mc
+    else:
+        if any(mags < mc):
+            raise ValueError("There are earthquakes below their respective "
+                             "completeness magnitude")
+
+    if n_m < min_num:
+        raise ValueError("n_m cannot be smaller than min_num")
+
+    if not isinstance(mags, np.ndarray):
+        raise ValueError("mags must be an array")
+    if not isinstance(times, np.ndarray):
+        raise ValueError("times must be an array")
+    if len(mags) != len(times):
+        raise ValueError("mags and times must have the same length")
+
+    if x_variable is None:
+        x_variable = np.arange(len(mags))
+    elif len(x_variable) != len(mags):
+        raise ValueError(
+            "x_variable must have the same length as magnitudes")
+    else:
+        idx_sort = np.argsort(x_variable)
+        mags = mags[idx_sort]
+        times = times[idx_sort]
+        mc = mc[idx_sort]
+        x_variable = x_variable[idx_sort]
+
+    if ax is None:
+        _, ax = plt.subplots()
+
+    if plot_technique == 'left':
+        idx_start = 0
+    elif plot_technique == 'midpoint':
+        idx_start = n_m // 2
+    elif plot_technique == 'right':
+        idx_start = n_m - 1
+    else:
+        raise ValueError(
+            "plot_technique must be 'left', 'midpoint', or 'right'")
+
+    # estimation
+    estimator = b_method()
+    b_values = np.ones(len(mags)) * np.nan
+    std_bs = np.ones(len(mags)) * np.nan
+    for ii in range(len(mags) - n_m + 1):
+        # runnning window of n_m magnitudes
+        mags_window = mags[ii:ii + n_m]
+        times_window = times[ii:ii + n_m]
+
+        # sort the magnitudes and times
+        idx = np.argsort(times_window)
+        mags_window = mags_window[idx]
+        times_window = times_window[idx]
+
+        # estimate the b-value
+        estimator.calculate(mags_window, mc=mc[ii], delta_m=delta_m, **kwargs)
+        if estimator.n < min_num:
+            b_values[idx_start + ii] = np.nan
+            std_bs[idx_start + ii] = np.nan
+        else:
+            b_values[idx_start + ii] = estimator.b_value
+            std_bs[idx_start + ii] = estimator.std
+
+    # plotting
+    ax.plot(x_variable, b_values, color=color, label=label)
+    error_factor = norm.ppf((1 + confidence) / 2)
+    ax.fill_between(
+        x_variable,
+        b_values - error_factor * std_bs,
+        b_values + error_factor * std_bs,
+        alpha=0.2,
+        color=color,
+        linewidth=0,
+    )
+
+    if label is not None:
+        ax.legend()
+    return ax