From 7eb62695be36a087c3f63e83f2f58244d87462bd Mon Sep 17 00:00:00 2001 From: Leila Mizrahi Date: Wed, 22 May 2024 14:12:01 +0200 Subject: [PATCH 1/7] docs: add automatic b-value calculation documentation --- docs/source/reference/analysis.md | 20 ++++++++++++++++++++ docs/source/reference/index.md | 1 + seismostats/analysis/__init__.py | 8 ++++++++ 3 files changed, 29 insertions(+) create mode 100644 docs/source/reference/analysis.md diff --git a/docs/source/reference/analysis.md b/docs/source/reference/analysis.md new file mode 100644 index 0000000..af00d75 --- /dev/null +++ b/docs/source/reference/analysis.md @@ -0,0 +1,20 @@ +# Analysis + +```{eval-rst} +.. currentmodule:: seismostats +``` + +## Estimating b-Values + +```{eval-rst} +.. autosummary:: + :toctree: api/ + + estimate_b + 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_2012 +``` \ 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..3894332 100644 --- a/seismostats/analysis/__init__.py +++ b/seismostats/analysis/__init__.py @@ -0,0 +1,8 @@ +# flake8: noqa + +from seismostats.analysis.estimate_beta import (estimate_b_kijko_smit_2012, + estimate_b_laplace, + estimate_b_positive, + estimate_b_tinti, + estimate_b_utsu, + estimate_b_weichert) From e9e515758387d609d162f635b21236ddf2c77316 Mon Sep 17 00:00:00 2001 From: Leila Mizrahi Date: Wed, 22 May 2024 14:19:28 +0200 Subject: [PATCH 2/7] doc: add mc methods and shi and bolt confidence --- docs/source/reference/analysis.md | 10 ++++++++++ seismostats/analysis/__init__.py | 1 + 2 files changed, 11 insertions(+) diff --git a/docs/source/reference/analysis.md b/docs/source/reference/analysis.md index af00d75..e05e8f9 100644 --- a/docs/source/reference/analysis.md +++ b/docs/source/reference/analysis.md @@ -11,10 +11,20 @@ :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_2012 +``` + +## 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/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py index 3894332..0ba67c5 100644 --- a/seismostats/analysis/__init__.py +++ b/seismostats/analysis/__init__.py @@ -6,3 +6,4 @@ estimate_b_tinti, estimate_b_utsu, estimate_b_weichert) +from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature From 751e5bb47a0f80f3310ecd49c87e4f4cf901cabc Mon Sep 17 00:00:00 2001 From: Leila Mizrahi Date: Wed, 22 May 2024 15:08:23 +0200 Subject: [PATCH 3/7] doc: formatting of the documentation, renaming kijko_smit_2012 --- docs/source/reference/analysis.md | 2 +- seismostats/analysis/__init__.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/reference/analysis.md b/docs/source/reference/analysis.md index e05e8f9..32a18bb 100644 --- a/docs/source/reference/analysis.md +++ b/docs/source/reference/analysis.md @@ -17,7 +17,7 @@ analysis.estimate_b_utsu analysis.estimate_b_laplace analysis.estimate_b_weichert - analysis.estimate_b_kijko_smit_2012 + analysis.estimate_b_kijko_smit ``` ## Estimating the Completeness Magnitude diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py index 0ba67c5..c8d6e2e 100644 --- a/seismostats/analysis/__init__.py +++ b/seismostats/analysis/__init__.py @@ -1,6 +1,6 @@ # flake8: noqa -from seismostats.analysis.estimate_beta import (estimate_b_kijko_smit_2012, +from seismostats.analysis.estimate_beta import (estimate_b_kijko_smit, estimate_b_laplace, estimate_b_positive, estimate_b_tinti, From c76bdf28f0ad5599fb1088f4467c7e40ab0116cf Mon Sep 17 00:00:00 2001 From: Leila Mizrahi Date: Wed, 22 May 2024 15:11:31 +0200 Subject: [PATCH 4/7] doc: formatting of the documentation, renaming kijko_smit_2012 --- seismostats/analysis/estimate_beta.py | 52 ++++++++++++++------------- seismostats/analysis/estimate_mc.py | 32 ++++++++++++----- 2 files changed, 52 insertions(+), 32 deletions(-) diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py index 58f8e3a..59195eb 100644 --- a/seismostats/analysis/estimate_beta.py +++ b/seismostats/analysis/estimate_beta.py @@ -33,14 +33,15 @@ 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) + + - '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 @@ -116,11 +117,12 @@ def estimate_b_tinti( b_parameter: str = "b_value", return_std: bool = False, ) -> float | tuple[float, float]: - """returns the maximum likelihood beta + """returns the maximum likelihood beta. + 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 @@ -173,6 +175,7 @@ def estimate_b_utsu( return_std: bool = False, ) -> float | tuple[float, float]: """returns the maximum likelihood beta + Source: Utsu 1965 (Geophysical bulletin of the Hokkaido University, vol 13, pp 99-103) @@ -294,7 +297,7 @@ def estimate_b_laplace( 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 + Magnitudes (this has a little less variance than the :func:`estimate_b_positive` method) Source: @@ -373,7 +376,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 +560,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, @@ -571,9 +574,10 @@ def estimate_b_kijko_smit_2012( 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 +589,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 +601,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 @@ -683,7 +687,7 @@ def shi_bolt_confidence( 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 + which parameter is given) according to shi and bolt 1982 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..d63187a 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) @@ -158,6 +158,13 @@ def mc_ks( Estimate the completeness magnitude (mc) 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 +285,17 @@ 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). + Estimate the completeness magnitude (mc) by maximum curvature. + + 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 From 4230cf2caa262685b225289890c2a93edffc8d00 Mon Sep 17 00:00:00 2001 From: Leila Mizrahi Date: Thu, 30 May 2024 09:37:12 +0200 Subject: [PATCH 5/7] documentation: homogenize the short description of each function to start with 'Return ...' --- seismostats/analysis/estimate_beta.py | 81 ++++++++++++++++----------- 1 file changed, 48 insertions(+), 33 deletions(-) diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py index 59195eb..533e700 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 @@ -34,17 +36,23 @@ def estimate_b( method: method to use for estimation of beta/b-value. Options are: - - 'tinti', 'utsu' (these are the classic estimators. 'tinti' + - '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) + 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) + 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, caution, can take long time to compute) + 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 + 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 @@ -117,11 +125,16 @@ 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 + - Tinti and Mulargia 1987 (Bulletin of the Seismological Society of America, 77(6), 2125-2134.) Args: @@ -174,7 +187,7 @@ 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 @@ -242,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 @@ -296,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 :func:`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 @@ -342,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 @@ -568,10 +584,9 @@ def estimate_b_kijko_smit( 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. and Smit, A., 2012. Extension of the Aki‐Utsu b‐value @@ -686,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 From 8a1e8b3c68c932110dbbad50ab0cdb4dbfaa804b Mon Sep 17 00:00:00 2001 From: Leila Mizrahi Date: Thu, 30 May 2024 09:41:50 +0200 Subject: [PATCH 6/7] documentation: start function descriptions with 'Return ...' --- seismostats/analysis/estimate_mc.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/seismostats/analysis/estimate_mc.py b/seismostats/analysis/estimate_mc.py index d63187a..0e4e417 100644 --- a/seismostats/analysis/estimate_mc.py +++ b/seismostats/analysis/estimate_mc.py @@ -155,7 +155,8 @@ 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: @@ -285,7 +286,8 @@ def mc_max_curvature( correction_factor: float = 0.2, ) -> float: """ - Estimate the completeness magnitude (mc) by maximum curvature. + Return the completeness magnitude (mc) estimate + using the maximum curvature method. Source: - Wiemer, S. and Wyss, M., 2000. Minimum magnitude of completeness From f62fba8465219586c776879a462e536d5acca618 Mon Sep 17 00:00:00 2001 From: Leila Mizrahi Date: Thu, 30 May 2024 09:51:28 +0200 Subject: [PATCH 7/7] fix: trailing whitespaces --- seismostats/analysis/estimate_beta.py | 4 ++-- seismostats/analysis/estimate_mc.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py index 533e700..7cc9896 100644 --- a/seismostats/analysis/estimate_beta.py +++ b/seismostats/analysis/estimate_beta.py @@ -360,7 +360,7 @@ def estimate_b_weichert( b_parameter: str = "b_value", ) -> tuple[float, float, float, float, float]: """Return the b-value estimate calculated using the - Weichert (1980) algorithm, for the case of unequal + Weichert (1980) algorithm, for the case of unequal completeness periods for different magnitude values. Source: @@ -585,7 +585,7 @@ def estimate_b_kijko_smit( b_parameter: str = 'b_value' ) -> tuple[float, float, float, float]: """Return the b-value estimate calculated using the - Kijko and Smit (2012) algorithm for the case of unequal + Kijko and Smit (2012) algorithm for the case of unequal completeness periods for different magnitude values. Source: diff --git a/seismostats/analysis/estimate_mc.py b/seismostats/analysis/estimate_mc.py index 0e4e417..aae6efc 100644 --- a/seismostats/analysis/estimate_mc.py +++ b/seismostats/analysis/estimate_mc.py @@ -162,8 +162,8 @@ def mc_ks( 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 + - 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: @@ -295,7 +295,7 @@ def mc_max_curvature( 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 + catalogues: Estimating the magnitude of completeness and its uncertainty. Bulletin of the Seismological Society of America, 95(2), pp.684-698. Args: