diff --git a/docs/source/reference/analysis.md b/docs/source/reference/analysis.md new file mode 100644 index 0000000..32a18bb --- /dev/null +++ b/docs/source/reference/analysis.md @@ -0,0 +1,30 @@ +# Analysis + +```{eval-rst} +.. currentmodule:: seismostats +``` + +## Estimating b-Values + +```{eval-rst} +.. autosummary:: + :toctree: api/ + + estimate_b + shi_bolt_confidence + analysis.estimate_b_tinti + analysis.estimate_b_positive + analysis.estimate_b_utsu + analysis.estimate_b_laplace + analysis.estimate_b_weichert + analysis.estimate_b_kijko_smit +``` + +## Estimating the Completeness Magnitude +```{eval-rst} +.. autosummary:: + :toctree: api/ + + analysis.mc_ks + analysis.mc_max_curvature +``` \ No newline at end of file diff --git a/docs/source/reference/index.md b/docs/source/reference/index.md index 83be9ac..ddc62f7 100644 --- a/docs/source/reference/index.md +++ b/docs/source/reference/index.md @@ -5,4 +5,5 @@ catalog config +analysis ``` \ No newline at end of file diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py index e69de29..c8d6e2e 100644 --- a/seismostats/analysis/__init__.py +++ b/seismostats/analysis/__init__.py @@ -0,0 +1,9 @@ +# flake8: noqa + +from seismostats.analysis.estimate_beta import (estimate_b_kijko_smit, + estimate_b_laplace, + estimate_b_positive, + estimate_b_tinti, + estimate_b_utsu, + estimate_b_weichert) +from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py index 58f8e3a..7cc9896 100644 --- a/seismostats/analysis/estimate_beta.py +++ b/seismostats/analysis/estimate_beta.py @@ -17,8 +17,10 @@ def estimate_b( method="tinti", return_n: bool = False, ) -> float | tuple[float, float] | tuple[float, float, float]: - """returns the maximum likelihood beta or b-value. Method depends on the - input parameter 'method'. + """Return the maximum likelihood estimator for the Gutenberg-Richter + b-value or beta, given an array of magnitudes + and a completeness magnitude. + Estimation method depends on the input parameter ``method``. Args: magnitudes: vector of magnitudes, unsorted, already cutoff (no @@ -33,17 +35,24 @@ def estimate_b( above) is returned method: method to use for estimation of beta/b-value. Options are: - - 'tinti', 'utsu' (these are the classic estimators. 'tinti' - is the recommended one, as it is more accurate. It is also - the default method) - - 'positive' (this is b-positive, which applies the 'tinti' - method to the positive differences. To achieve the effect - of reduced STAI, the magnitudes must be ordered in time) - - 'laplace' (this is using the distribution of all - differences, caution, can take long time to compute) - return_n: if True the number of events used for the estimation is - returned. This is only relevant for the 'positive' method + - 'tinti', 'utsu' (these are the classic estimators, see + :func:`seismostats.analysis.estimate_b_tinti`, + :func:`seismostats.analysis.estimate_b_utsu`. 'tinti' + is the recommended one, as it is more accurate. It is also + the default method.) + - 'positive' (this is b-positive, which applies the 'tinti' + method to the positive differences, see + :func:`seismostats.analysis.estimate_b_positive`. + To achieve the effect + of reduced STAI, the magnitudes must be ordered in time.) + - 'laplace' (this is using the distribution of all + differences, see + :func:`seismostats.analysis.estimate_b_laplace`. + Caution: it can take long time to compute.) + + return_n: If True, the number of events used for the estimation is + returned. This is only relevant for the 'positive' method. Returns: b: maximum likelihood beta or b-value, depending on value of @@ -116,11 +125,17 @@ def estimate_b_tinti( b_parameter: str = "b_value", return_std: bool = False, ) -> float | tuple[float, float]: - """returns the maximum likelihood beta + """Return the maximum likelihood b-value or beta for + an array of magnitudes and a completeness magnitude mc. + If the magnitudes are discretized, the discretization must be given in + ``delta_m``, so that the maximum likelihood estimator can be calculated + correctly. + + Source: - Aki 1965 (Bull. Earthquake research institute, vol 43, pp 237-239) - Tinti and Mulargia 1987 (Bulletin of the Seismological Society of - America, 77(6), 2125-2134.) + - Aki 1965 (Bull. Earthquake research institute, vol 43, pp 237-239) + - Tinti and Mulargia 1987 (Bulletin of the Seismological Society of + America, 77(6), 2125-2134.) Args: magnitudes: vector of magnitudes, unsorted, already cutoff (no @@ -172,7 +187,8 @@ def estimate_b_utsu( b_parameter: str = "b_value", return_std: bool = False, ) -> float | tuple[float, float]: - """returns the maximum likelihood beta + """Return the maximum likelihood b-value or beta. + Source: Utsu 1965 (Geophysical bulletin of the Hokkaido University, vol 13, pp 99-103) @@ -239,24 +255,26 @@ def estimate_b_positive( return_std: bool = False, return_n: bool = False, ) -> float | tuple[float, float] | tuple[float, float, float]: - """returns the b-value estimation using the positive differences of the - Magnitudes + """Return the b-value estimate calculated using the + positive differences between consecutive magnitudes. Source: Van der Elst 2021 (J Geophysical Research: Solid Earth, Vol 126, Issue 2) Args: - magnitudes: vector of magnitudes differences, sorted in time (first - entry is the earliest earthquake) - delta_m: discretization of magnitudes. default is no discretization + magnitudes: array of magnitudes, sorted in time (first + entry is the earliest earthquake). + To achieve the effect + of reduced STAI, the magnitudes must be ordered in time. + delta_m: discretization of magnitudes. default is no discretization. b_parameter:either 'b-value', then the corresponding value of the Gutenberg-Richter law is returned, otherwise 'beta' from the - exponential distribution [p(M) = exp(-beta*(M-mc))] + exponential distribution [p(M) = exp(-beta*(M-mc))]. return_std: if True the standard deviation of beta/b-value (see above) - is returned + is returned. return_n: if True the number of events used for the estimation is - returned + returned. Returns: b: maximum likelihood beta or b-value, depending on value of @@ -293,9 +311,11 @@ def estimate_b_laplace( return_std: bool = False, return_n: bool = False, ) -> float | tuple[float, float]: - """returns the b-value estimation using the all the differences of the - Magnitudes (this has a little less variance than the estimate_b_positive - method) + """Return the b-value estimate calculated using + all the differences between magnitudes. + (This has a little less variance than the + :func:`seismostats.analysis.estimate_b_positive` + method.) Source: Van der Elst 2021 (J Geophysical Research: Solid Earth, Vol 126, Issue @@ -339,10 +359,9 @@ def estimate_b_weichert( delta_m: float = 0.1, b_parameter: str = "b_value", ) -> tuple[float, float, float, float, float]: - """applies the Weichert (1980) algorithm for estimation of the - Gutenberg-Richter magnitude-frequency distribution parameters in - the case of unequal completeness periods for different magnitude - values. + """Return the b-value estimate calculated using the + Weichert (1980) algorithm, for the case of unequal + completeness periods for different magnitude values. Source: Weichert (1980), Estimation of the earthquake recurrence parameters @@ -373,7 +392,7 @@ def estimate_b_weichert( Gutenberg-Richter law is returned, otherwise 'beta' from the exponential distribution [p(M) = exp(-beta*(M-mc))] - Returns:( + Returns: b_parameter: maximum likelihood point estimate of 'b-value' or 'beta' std_b_parameter: standard error of b_parameter rate_at_lmc: maximum likelihood point estimate of earthquake rate @@ -557,7 +576,7 @@ def _weichert_objective_function( return np.abs(left - right) -def estimate_b_kijko_smit_2012( +def estimate_b_kijko_smit( magnitudes: np.ndarray, dates: list[np.datetime64], completeness_table: np.ndarray, @@ -565,15 +584,15 @@ def estimate_b_kijko_smit_2012( delta_m: float = 0.1, b_parameter: str = 'b_value' ) -> tuple[float, float, float, float]: - """ applies the Kijko and Smit (2012) algorithm for estimation of the - Gutenberg-Richter magnitude-frequency distribution parameters in - the case of unequal completeness periods for different magnitude - values. + """Return the b-value estimate calculated using the + Kijko and Smit (2012) algorithm for the case of unequal + completeness periods for different magnitude values. Source: - Kijko A., Smit A. (2012), Extension of the Aki‐Utsu b‐Value Estimator - for Incomplete Catalogs, Bulletin of the Seismological Society of - America Vol 102, No. 3, pp. 1283–1287 + Kijko, A. and Smit, A., 2012. Extension of the Aki‐Utsu b‐value + estimator for incomplete catalogs. Bulletin of the Seismological + Society of America, 102(3), pp.1283-1287. + Args: magnitudes: vector of earthquake magnitudes @@ -585,10 +604,10 @@ def estimate_b_kijko_smit_2012( the first column are considered detected. An example is given below: - np.array([[ 3.95, 1980], - [ 4.95, 1920], - [ 5.95, 1810], - [ 6.95, 1520]]) + >>> np.array([[ 3.95, 1980], + ... [ 4.95, 1920], + ... [ 5.95, 1810], + ... [ 6.95, 1520]]) last_year: last year of observation (the default is None, in which case it is set to the latest year in years). @@ -597,7 +616,7 @@ def estimate_b_kijko_smit_2012( Gutenberg-Richter law is returned, otherwise 'beta' from the exponential distribution [p(M) = exp(-beta*(M-mc))] - Returns:( + Returns: b_parameter: maximum likelihood point estimate of 'b-value' or 'beta' std_b_parameter: standard error of b_parameter rate_at_lmc: maximum likelihood point estimate of earthquake rate @@ -682,8 +701,8 @@ def shi_bolt_confidence( b: float | None = None, b_parameter: str = "b_value", ) -> float: - """calculates the confidence limit of the b_value or beta (depending on - which parameter is given) according to shi and bolt 1982 + """Return the Shi and Bolt (1982) confidence limit of the b-value or + beta. Source: Shi and Bolt, BSSA, Vol. 72, No. 5, pp. 1677-1687, October 1982 diff --git a/seismostats/analysis/estimate_mc.py b/seismostats/analysis/estimate_mc.py index 36c4535..aae6efc 100644 --- a/seismostats/analysis/estimate_mc.py +++ b/seismostats/analysis/estimate_mc.py @@ -2,15 +2,15 @@ for the estimation of the completeness magnitude. """ +import warnings + import numpy as np import pandas as pd -import warnings from seismostats.analysis.estimate_beta import estimate_b -from seismostats.utils.binning import get_fmd, bin_to_precision -from seismostats.utils.simulate_distributions import ( - simulated_magnitudes_binned, -) +from seismostats.utils.binning import bin_to_precision, get_fmd +from seismostats.utils.simulate_distributions import \ + simulated_magnitudes_binned def cdf_discrete_GR( @@ -127,7 +127,7 @@ def ks_test_gr( ) for ii in range(n): - simulated = simulated_all[n_sample * ii : n_sample * (ii + 1)] + simulated = simulated_all[n_sample * ii: n_sample * (ii + 1)] _, y_th = cdf_discrete_GR(simulated, mc=mc, delta_m=delta_m, beta=beta) _, y_emp = empirical_cdf(simulated) @@ -155,9 +155,17 @@ def mc_ks( n: int = 10000, ) -> tuple[np.ndarray, list[float], np.ndarray, float | None, float | None]: """ - Estimate the completeness magnitude (mc) for a given list of completeness + Return the completeness magnitude (mc) estimate + for a given list of completeness magnitudes using the K-S distance method. + Source: + - Clauset, A., Shalizi, C.R. and Newman, M.E., 2009. Power-law + distributions in empirical data. SIAM review, 51(4), pp.661-703. + - Mizrahi, L., Nandan, S. and Wiemer, S., 2021. The effect of + declustering on the size distribution of mainshocks. Seismological + Society of America, 92(4), pp.2333-2342. + Args: sample: Magnitudes to test delta_m: Magnitude bins (sample has to be rounded to bins @@ -278,8 +286,18 @@ def mc_max_curvature( correction_factor: float = 0.2, ) -> float: """ - Estimate the completeness magnitude (mc) by maximum curvature (Wiemer and - Wyss 2000, Woessner and Wiemer 2005). + Return the completeness magnitude (mc) estimate + using the maximum curvature method. + + Source: + - Wiemer, S. and Wyss, M., 2000. Minimum magnitude of completeness + in earthquake catalogs: Examples from Alaska, the western United + States, and Japan. Bulletin of the Seismological Society of America, + 90(4), pp.859-869. + - Woessner, J. and Wiemer, S., 2005. Assessing the quality of earthquake + catalogues: Estimating the magnitude of completeness and its + uncertainty. + Bulletin of the Seismological Society of America, 95(2), pp.684-698. Args: sample: Magnitudes to test delta_m: Magnitude bins (sample has to be rounded to bins