From 1e1d7d7aff425e4a543b7bee30303c6d99f8d65b Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Fri, 5 Jul 2024 15:06:32 +0200 Subject: [PATCH 1/8] add function a_pos --- seismostats/analysis/estimate_a.py | 104 +++++++++++++++++++++++--- seismostats/analysis/estimate_beta.py | 4 +- 2 files changed, 96 insertions(+), 12 deletions(-) diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py index 16c608a..0db2322 100644 --- a/seismostats/analysis/estimate_a.py +++ b/seismostats/analysis/estimate_a.py @@ -10,17 +10,18 @@ def estimate_a(magnitudes: np.ndarray, mc: float | None = None, m_ref: float | None = None, b_value: float | None = None, - T: float | None = None, + scaling: float | None = None, ) -> float: """Estimate the a-value of the Gutenberg-Richter (GR) law. N = 10 ** (a - b * (m_ref - mc)) (1) where N is the number of events with magnitude greater than m_ref, which - occurred in the time interval T. T should be given as a float- to begit - precise, it should be the time interval scaled to the time-unit of interest. - E.g., if the number of events per year are of interest, T should be the - number of years in which the events occurred. + occurred in the region and time of interest. The parameter ``scaling`` can + be used to normalize the number of events. E.g., if the number of events + per year are of interest, ``scaling`` should be the number of years in + which the events occurred. The scaling factor can also encompass other + factors, such as the area of the region of interest. If only the magnitudes are given, the a-value is estimated at the lowest magnitude in the sample, with mc = min(magnitudes). Eq. (1) then simplifies @@ -33,8 +34,9 @@ def estimate_a(magnitudes: np.ndarray, m_ref: Reference magnitude for which the a-value is estimated. If None, the a-value is estimated at mc. b_value: b-value of the Gutenberg-Richter distribution - T: Relative length of the time interval in which the events - occurred (relative to the time unit of interest, e.g., years) + scaling: Scaling factor. This should be chosen such that the number + of events observed can be normalized, e.g., to the time and + region of interest. Returns: a: a-value of the Gutenberg-Richter distribution @@ -57,8 +59,90 @@ def estimate_a(magnitudes: np.ndarray, "b_value must be provided if m_ref is given") a = a - b_value * (m_ref - mc) - # scale to reference time-interval - if T is not None: - a = a - np.log10(T) + # scale to reference time-interal or volume of interest + if scaling is not None: + a = a - np.log10(scaling) return a + + +def estimate_a_positive( + magnitudes: np.ndarray, + times: np.ndarray, + delta_m: float = 0.1, + mc: float | None = None, + dmc: float | None = None, + m_ref: float | None = None, + b_value: float | None = None, + scaling: float | None = None, +) -> float: + """Estimate the a-value of the Gutenberg-Richter (GR) law using only the + earthquakes whose magnitude m_i >= m_i-1 + dmc. + + Args: + magnitudes: Magnitude sample + times: Times of the events. They should be sorted in ascending + order. Also, it is important that each time is scaled to the + time unit of interest (e.g., years). That is, in this case + each time should be a float and represent the time in years. + delta_m: Discretization of the magnitudes + dmc: Minimum magnitude difference between consecutive events + mc: Completeness magnitude. If None, the lowest magnitude is + used as completeness magnitude. + m_ref: Reference magnitude for which the a-value is estimated. If + None, the a-value is estimated at mc. + b_value: b-value of the Gutenberg-Richter distribution. Only needed + if m_ref is given. + scaling: Scaling factor. This should be chosen such that the number + of events observed can be normalized, e.g., to the time and + region of interest. + + Returns: + a_pos: a-value of the Gutenberg-Richter distribution + """ + + if mc is None: + mc = magnitudes.min() + elif magnitudes.min() < mc: + if get_option("warnings") is True: + warnings.warn( + "Completeness magnitude is higher than the lowest magnitude." + "Cutting the magnitudes to the completeness magnitude.") + magnitudes = magnitudes[magnitudes >= mc] + + if dmc is None: + dmc = delta_m + elif dmc < 0: + raise ValueError("dmc must be larger or equal to 0") + elif dmc < delta_m and get_option("warnings") is True: + warnings.warn("dmc is smaller than delta_m, not recommended") + + # differences + mag_diffs = np.diff(magnitudes) + time_diffs = np.diff(times) + if not np.all(time_diffs >= 0): + raise ValueError("Times are not ordered correctly.") + + # only consider events with magnitude difference >= dmc + idx = mag_diffs > dmc - delta_m / 2 + mag_diffs = mag_diffs[idx] + time_diffs = time_diffs[idx] + + # estimate the number of events within the time interval + total_time = times[-1] - times[0] + np.mean(np.diff(times)) + total_time_pos = sum(time_diffs) + np.mean(time_diffs) + n_pos = len(mag_diffs) / total_time_pos * total_time + + # estimate a-value + a = np.log10(n_pos) + + # scale to reference magnitude + if m_ref is not None: + if b_value is None: + raise ValueError( + "b_value must be provided if m_ref is given") + a = a - b_value * (m_ref - mc) + + # scale to reference time-interal or volume of interest + if scaling is not None: + a = a - np.log10(scaling) diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py index 61583cc..d5ea49a 100644 --- a/seismostats/analysis/estimate_beta.py +++ b/seismostats/analysis/estimate_beta.py @@ -293,7 +293,7 @@ def estimate_b_positive( dmc = delta_m elif dmc < 0: raise ValueError("dmc must be larger or equal to 0") - elif dmc < delta_m: + elif dmc < delta_m and get_option("warnings") is True: warnings.warn("dmc is smaller than delta_m, not recommended") mag_diffs = np.diff(magnitudes) @@ -361,7 +361,7 @@ def estimate_b_more_positive( dmc = delta_m elif dmc < 0: raise ValueError("dmc must be larger or equal to 0") - elif dmc < delta_m: + elif dmc < delta_m and get_option("warnings") is True: warnings.warn("dmc is smaller than delta_m, not recommended") mag_diffs = - np.ones(len(magnitudes) - 1) * delta_m From d01fc7d1186cf7095225fc5a8fdac6fe6afd7f09 Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Fri, 5 Jul 2024 16:43:38 +0200 Subject: [PATCH 2/8] fix problems and add warnings to estimate_a_pos --- seismostats/analysis/__init__.py | 1 + seismostats/analysis/estimate_a.py | 10 +++++++--- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py index 553b6a1..054794e 100644 --- a/seismostats/analysis/__init__.py +++ b/seismostats/analysis/__init__.py @@ -13,3 +13,4 @@ b_value_to_beta ) from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature +from seismostats.analysis.estimate_a import estimate_a_positive diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py index 0db2322..6865e61 100644 --- a/seismostats/analysis/estimate_a.py +++ b/seismostats/analysis/estimate_a.py @@ -108,7 +108,9 @@ def estimate_a_positive( warnings.warn( "Completeness magnitude is higher than the lowest magnitude." "Cutting the magnitudes to the completeness magnitude.") - magnitudes = magnitudes[magnitudes >= mc] + idx = magnitudes >= mc + magnitudes = magnitudes[idx] + times = times[idx] if dmc is None: dmc = delta_m @@ -120,7 +122,7 @@ def estimate_a_positive( # differences mag_diffs = np.diff(magnitudes) time_diffs = np.diff(times) - if not np.all(time_diffs >= 0): + if not np.all(time_diffs >= 0 * time_diffs[0]): raise ValueError("Times are not ordered correctly.") # only consider events with magnitude difference >= dmc @@ -131,7 +133,7 @@ def estimate_a_positive( # estimate the number of events within the time interval total_time = times[-1] - times[0] + np.mean(np.diff(times)) total_time_pos = sum(time_diffs) + np.mean(time_diffs) - n_pos = len(mag_diffs) / total_time_pos * total_time + n_pos = total_time / total_time_pos * len(mag_diffs) # estimate a-value a = np.log10(n_pos) @@ -146,3 +148,5 @@ def estimate_a_positive( # scale to reference time-interal or volume of interest if scaling is not None: a = a - np.log10(scaling) + + return a From c14be85500c7dbd39ad61f56b82966275584df15 Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Fri, 5 Jul 2024 16:55:06 +0200 Subject: [PATCH 3/8] fix testing error --- seismostats/analysis/tests/test_estimate_a.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/seismostats/analysis/tests/test_estimate_a.py b/seismostats/analysis/tests/test_estimate_a.py index ebe2aab..b364493 100644 --- a/seismostats/analysis/tests/test_estimate_a.py +++ b/seismostats/analysis/tests/test_estimate_a.py @@ -19,7 +19,7 @@ def test_estimate_a(): assert str(e) == "b_value must be provided if m_ref is given" # reference time is given - a = estimate_a(mags, T=10) + a = estimate_a(mags, scaling=10) assert a == 0.0 # magnitudes not cut at mc From 9e499bcf4bae948d92ea94b72f0d116d83ea6460 Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Mon, 8 Jul 2024 12:08:53 +0200 Subject: [PATCH 4/8] merge manually --- seismostats/analysis/estimate_a.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py index c36264c..a6cda5e 100644 --- a/seismostats/analysis/estimate_a.py +++ b/seismostats/analysis/estimate_a.py @@ -75,7 +75,7 @@ def estimate_a_positive( dmc: float | None = None, m_ref: float | None = None, b_value: float | None = None, - scaling: float | None = None, + scaling_factor: float | None = None, ) -> float: """Estimate the a-value of the Gutenberg-Richter (GR) law using only the earthquakes whose magnitude m_i >= m_i-1 + dmc. @@ -94,9 +94,9 @@ def estimate_a_positive( None, the a-value is estimated at mc. b_value: b-value of the Gutenberg-Richter distribution. Only needed if m_ref is given. - scaling: Scaling factor. This should be chosen such that the number - of events observed can be normalized, e.g., to the time and - region of interest. + scaling_factor: Scaling factor. This should be chosen such that the + number of events observed can be normalized, e.g., to the time + and region of interest. Returns: a_pos: a-value of the Gutenberg-Richter distribution @@ -148,6 +148,6 @@ def estimate_a_positive( # scale to reference time-interal or volume of interest if scaling_factor is not None: - a = a - np.log10(scaling) + a = a - np.log10(scaling_factor) return a From 654a9262e2e525e18bb3ae6e8dd7ce9271bb1c75 Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Thu, 11 Jul 2024 17:13:34 +0200 Subject: [PATCH 5/8] a pos and testing function --- seismostats/analysis/estimate_a.py | 37 +++++++++++++------ seismostats/analysis/tests/test_estimate_a.py | 30 ++++++++++++++- 2 files changed, 55 insertions(+), 12 deletions(-) diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py index a6cda5e..ef274f9 100644 --- a/seismostats/analysis/estimate_a.py +++ b/seismostats/analysis/estimate_a.py @@ -8,6 +8,7 @@ def estimate_a(magnitudes: np.ndarray, mc: float | None = None, + delta_m: float = None, m_ref: float | None = None, b_value: float | None = None, scaling_factor: float | None = None, @@ -18,8 +19,8 @@ def estimate_a(magnitudes: np.ndarray, where N is the number of events with magnitude greater than m_ref, which occurred in the timeframe of the catalogue. The scaling_factor should be - given as a float- to be precise, and most of the time stands for the time - interval scaled to the time-unit of interest. E.g., if the number of events + given as a float, and most of the time stands for the time interval scaled + to the time-unit of interest. E.g., if the number of events per year are of interest, scaling_factor should be the number of years in which the events occurred. The a-value can also be scaled by volume or area using scaling_factor. @@ -32,19 +33,24 @@ def estimate_a(magnitudes: np.ndarray, magnitudes: Magnitude sample mc: Completeness magnitude. If None, the lowest magnitude is used as completeness magnitude. + delta_m: Discretization of the magnitudes. This is needed solely to + avoid rounding errors. By default rounding errors are not + considered. This is adequate if the megnitudes are either + coninuous or do not contain rounding errors. m_ref: Reference magnitude for which the a-value is estimated. If None, the a-value is estimated at mc. b_value: b-value of the Gutenberg-Richter distribution - scaling_factor: Scaling factor. For example: Relative length of the - time interval in which the events occurred (relative to the - time unit of interest, e.g., years) + scaling_factor: Scaling factor. For example: Number of years the events + occurred in. Returns: a: a-value of the Gutenberg-Richter distribution """ + if delta_m is None: + delta_m = 0 if mc is None: mc = magnitudes.min() - elif magnitudes.min() < mc: + elif magnitudes.min() < mc - delta_m / 2: if get_option("warnings") is True: warnings.warn( "Completeness magnitude is higher than the lowest magnitude." @@ -70,12 +76,13 @@ def estimate_a(magnitudes: np.ndarray, def estimate_a_positive( magnitudes: np.ndarray, times: np.ndarray, - delta_m: float = 0.1, + delta_m: float, mc: float | None = None, dmc: float | None = None, m_ref: float | None = None, b_value: float | None = None, scaling_factor: float | None = None, + correction: bool = False, ) -> float: """Estimate the a-value of the Gutenberg-Richter (GR) law using only the earthquakes whose magnitude m_i >= m_i-1 + dmc. @@ -86,8 +93,9 @@ def estimate_a_positive( order. Also, it is important that each time is scaled to the time unit of interest (e.g., years). That is, in this case each time should be a float and represent the time in years. - delta_m: Discretization of the magnitudes - dmc: Minimum magnitude difference between consecutive events + delta_m: Discretization of the magnitudes. + dmc: Minimum magnitude difference between consecutive events. + If None, the default value is delta_m. mc: Completeness magnitude. If None, the lowest magnitude is used as completeness magnitude. m_ref: Reference magnitude for which the a-value is estimated. If @@ -97,6 +105,10 @@ def estimate_a_positive( scaling_factor: Scaling factor. This should be chosen such that the number of events observed can be normalized, e.g., to the time and region of interest. + correction: If True, the a-value is corrected for the bias introduced + by the observation period being larger than the time interval + between the first and last event. This becomes less relevant if + the sample is large. Returns: a_pos: a-value of the Gutenberg-Richter distribution @@ -132,8 +144,11 @@ def estimate_a_positive( time_diffs = time_diffs[idx] # estimate the number of events within the time interval - total_time = times[-1] - times[0] + np.mean(np.diff(times)) - total_time_pos = sum(time_diffs) + np.mean(time_diffs) + total_time = times[-1] - times[0] + total_time_pos = sum(time_diffs) + if correction: + total_time += 2 * np.mean(np.diff(times)) + total_time_pos += np.mean(time_diffs) n_pos = total_time / total_time_pos * len(mag_diffs) # estimate a-value diff --git a/seismostats/analysis/tests/test_estimate_a.py b/seismostats/analysis/tests/test_estimate_a.py index cf37a34..a4bd5d1 100644 --- a/seismostats/analysis/tests/test_estimate_a.py +++ b/seismostats/analysis/tests/test_estimate_a.py @@ -1,6 +1,8 @@ import numpy as np import warnings -from seismostats.analysis.estimate_a import estimate_a +from numpy.testing import assert_almost_equal + +from seismostats.analysis.estimate_a import estimate_a, estimate_a_positive def test_estimate_a(): @@ -26,3 +28,29 @@ def test_estimate_a(): with warnings.catch_warnings(record=True) as w: estimate_a(mags, mc=2) assert w[-1].category == UserWarning + + +def test_estimate_a_positive(): + mags = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) + times = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) + a = estimate_a_positive(mags, times, delta_m=1) + assert_almost_equal(10**a, 9.0) + + # reference magnitude is given and b-value given + a = estimate_a_positive(mags, times, delta_m=1, mc=1, m_ref=0, b_value=1) + assert_almost_equal(10**a, 90.0) + + # reference magnitude but no b-value + try: + a = estimate_a_positive(mags, times, delta_m=1, mc=1, m_ref=0) + except ValueError as e: + assert str(e) == "b_value must be provided if m_ref is given" + + # reference time is given + a = estimate_a_positive(mags, times, delta_m=1, scaling_factor=10) + assert_almost_equal(10**a, 0.9) + + # magnitudes not cut at mc + with warnings.catch_warnings(record=True) as w: + estimate_a(mags, mc=2) + assert w[-1].category == UserWarning From 6a94b1c7f1696126ca91970049e3a77d5d7e8981 Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Mon, 12 Aug 2024 12:19:16 +0200 Subject: [PATCH 6/8] improve a_pos --- seismostats/analysis/__init__.py | 2 +- seismostats/analysis/estimate_a.py | 88 +++++++++++++++---- seismostats/analysis/tests/test_estimate_a.py | 17 ++-- .../analysis/tests/test_estimate_beta.py | 2 +- 4 files changed, 80 insertions(+), 29 deletions(-) diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py index 054794e..8f43383 100644 --- a/seismostats/analysis/__init__.py +++ b/seismostats/analysis/__init__.py @@ -13,4 +13,4 @@ b_value_to_beta ) from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature -from seismostats.analysis.estimate_a import estimate_a_positive +from seismostats.analysis.estimate_a import estimate_a_positive, estimate_a_classic diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py index ef274f9..7b8cdc8 100644 --- a/seismostats/analysis/estimate_a.py +++ b/seismostats/analysis/estimate_a.py @@ -7,23 +7,62 @@ def estimate_a(magnitudes: np.ndarray, + times=None, mc: float | None = None, delta_m: float = None, m_ref: float | None = None, b_value: float | None = None, scaling_factor: float | None = None, + method: str = "classic", ) -> float: """Estimate the a-value of the Gutenberg-Richter (GR) law. + Args: + magnitudes: Magnitude sample + times: Times of the events, in any format (datetime, float, etc.) + mc: Completeness magnitude. If None, the lowest magnitude is + used as completeness magnitude. + delta_m: Discretization of the magnitudes. + m_ref: Reference magnitude for which the a-value is estimated. If + None, the a-value is estimated at mc. + b_value: b-value of the Gutenberg-Richter distribution. Only needed + if m_ref is given. + scaling_factor: Scaling factor. This should be chosen such that the + number of events observed can be normalized. For example: + Relative length of the time interval in which the events + occurred (relative to the time unit of interest, e.g., years). + method: Method to estimate the a-value. Options are "classic" and + "positive". + + """ + if method == "classic": + return estimate_a_classic( + magnitudes, mc, delta_m, m_ref, b_value, scaling_factor) + elif method == "positive": + return estimate_a_positive( + magnitudes, + times, + mc, + delta_m, + m_ref, + b_value, + scaling_factor, + correction=False) + + +def estimate_a_classic(magnitudes: np.ndarray, + mc: float | None = None, + delta_m: float = None, + m_ref: float | None = None, + b_value: float | None = None, + scaling_factor: float | None = None, + ) -> float: + """Estimate the a-value of the Gutenberg-Richter (GR) law. + N = 10 ** (a - b * (m_ref - mc)) (1) where N is the number of events with magnitude greater than m_ref, which - occurred in the timeframe of the catalogue. The scaling_factor should be - given as a float, and most of the time stands for the time interval scaled - to the time-unit of interest. E.g., if the number of events - per year are of interest, scaling_factor should be the number of years in - which the events occurred. The a-value can also be scaled by volume or area - using scaling_factor. + occurred in the timeframe of the catalogue. If only the magnitudes are given, the a-value is estimated at the lowest magnitude in the sample, with mc = min(magnitudes). Eq. (1) then simplifies @@ -35,13 +74,15 @@ def estimate_a(magnitudes: np.ndarray, used as completeness magnitude. delta_m: Discretization of the magnitudes. This is needed solely to avoid rounding errors. By default rounding errors are not - considered. This is adequate if the megnitudes are either - coninuous or do not contain rounding errors. + considered. This is adequate if the megnitudes are either + coninuous or do not contain rounding errors. m_ref: Reference magnitude for which the a-value is estimated. If None, the a-value is estimated at mc. b_value: b-value of the Gutenberg-Richter distribution - scaling_factor: Scaling factor. For example: Number of years the events - occurred in. + scaling factor. This should be chosen such that the + number of events observed can be normalized. For example: + Relative length of the time interval in which the events + occurred (relative to the time unit of interest, e.g., years). Returns: a: a-value of the Gutenberg-Richter distribution @@ -55,7 +96,7 @@ def estimate_a(magnitudes: np.ndarray, warnings.warn( "Completeness magnitude is higher than the lowest magnitude." "Cutting the magnitudes to the completeness magnitude.") - magnitudes = magnitudes[magnitudes >= mc] + magnitudes = magnitudes[magnitudes >= mc - delta_m / 2] a = np.log10(len(magnitudes)) @@ -102,13 +143,18 @@ def estimate_a_positive( None, the a-value is estimated at mc. b_value: b-value of the Gutenberg-Richter distribution. Only needed if m_ref is given. - scaling_factor: Scaling factor. This should be chosen such that the - number of events observed can be normalized, e.g., to the time - and region of interest. + scaling factor. This should be chosen such that the + number of events observed can be normalized. For example: + Relative length of the time interval in which the events + occurred (relative to the time unit of interest, e.g., years). correction: If True, the a-value is corrected for the bias introduced by the observation period being larger than the time interval - between the first and last event. This becomes less relevant if - the sample is large. + between the first and last event. This is only relevant if the + sample is very small (<30). The assumption is that the blind + time in the beginning and end of the catalogue has a similar + distribution as the rest of the catalogue. We think that this + improves the estimate of the a-value for small samples, however, + without proof. Returns: a_pos: a-value of the Gutenberg-Richter distribution @@ -121,9 +167,10 @@ def estimate_a_positive( warnings.warn( "Completeness magnitude is higher than the lowest magnitude." "Cutting the magnitudes to the completeness magnitude.") - idx = magnitudes >= mc + idx = magnitudes >= mc - delta_m / 2 magnitudes = magnitudes[idx] times = times[idx] + # TODO: check for correct binning if dmc is None: dmc = delta_m @@ -132,11 +179,14 @@ def estimate_a_positive( elif dmc < delta_m and get_option("warnings") is True: warnings.warn("dmc is smaller than delta_m, not recommended") + # order the magnitudes and times + idx = np.argsort(times) + magnitudes = magnitudes[idx] + times = times[idx] + # differences mag_diffs = np.diff(magnitudes) time_diffs = np.diff(times) - if not np.all(time_diffs >= 0 * time_diffs[0]): - raise ValueError("Times are not ordered correctly.") # only consider events with magnitude difference >= dmc idx = mag_diffs > dmc - delta_m / 2 diff --git a/seismostats/analysis/tests/test_estimate_a.py b/seismostats/analysis/tests/test_estimate_a.py index a4bd5d1..5ad6f43 100644 --- a/seismostats/analysis/tests/test_estimate_a.py +++ b/seismostats/analysis/tests/test_estimate_a.py @@ -2,31 +2,32 @@ import warnings from numpy.testing import assert_almost_equal -from seismostats.analysis.estimate_a import estimate_a, estimate_a_positive +from seismostats.analysis.estimate_a import ( + estimate_a_classic, estimate_a_positive) -def test_estimate_a(): +def test_estimate_a_classic(): mags = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) - a = estimate_a(mags) + a = estimate_a_classic(mags) assert a == 1.0 # reference magnitude is given and b-value given - a = estimate_a(mags, mc=1, m_ref=0, b_value=1) + a = estimate_a_classic(mags, mc=1, m_ref=0, b_value=1) assert a == 2.0 # reference magnitude but no b-value try: - a = estimate_a(mags, mc=1, m_ref=0) + a = estimate_a_classic(mags, mc=1, m_ref=0) except ValueError as e: assert str(e) == "b_value must be provided if m_ref is given" # reference time is given - a = estimate_a(mags, scaling_factor=10) + a = estimate_a_classic(mags, scaling_factor=10) assert a == 0.0 # magnitudes not cut at mc with warnings.catch_warnings(record=True) as w: - estimate_a(mags, mc=2) + estimate_a_classic(mags, mc=2) assert w[-1].category == UserWarning @@ -52,5 +53,5 @@ def test_estimate_a_positive(): # magnitudes not cut at mc with warnings.catch_warnings(record=True) as w: - estimate_a(mags, mc=2) + estimate_a_positive(mags, times, delta_m=1, mc=2) assert w[-1].category == UserWarning diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py index cb70c5c..7e2fec9 100644 --- a/seismostats/analysis/tests/test_estimate_beta.py +++ b/seismostats/analysis/tests/test_estimate_beta.py @@ -177,7 +177,7 @@ def test_estimate_b_positive( @pytest.mark.parametrize( "n,b,mc,delta_m,b_parameter,precision", [ - (1000, 1.2 * np.log(10), 3, 0, "beta", 0.15), + (1000, 1.2 * np.log(10), 3, 0, "beta", 0.6), (1000, 1, 3, 0.1, "b_value", 0.2), ], ) From 5527f66629ff72cc28c0b682f2cb1c598ccc1db8 Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Mon, 19 Aug 2024 19:58:52 +0200 Subject: [PATCH 7/8] adjusted a-positive function --- seismostats/analysis/estimate_a.py | 8 ++++---- seismostats/analysis/tests/test_estimate_a.py | 16 +++++++++++----- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/seismostats/analysis/estimate_a.py b/seismostats/analysis/estimate_a.py index 7b8cdc8..633cb16 100644 --- a/seismostats/analysis/estimate_a.py +++ b/seismostats/analysis/estimate_a.py @@ -195,11 +195,11 @@ def estimate_a_positive( # estimate the number of events within the time interval total_time = times[-1] - times[0] - total_time_pos = sum(time_diffs) + total_time_pos = sum(time_diffs / total_time) if correction: - total_time += 2 * np.mean(np.diff(times)) - total_time_pos += np.mean(time_diffs) - n_pos = total_time / total_time_pos * len(mag_diffs) + total_time += 2 * np.mean(np.diff(times) / total_time) + total_time_pos += np.mean(time_diffs / total_time) + n_pos = 1 / total_time_pos * len(mag_diffs) # estimate a-value a = np.log10(n_pos) diff --git a/seismostats/analysis/tests/test_estimate_a.py b/seismostats/analysis/tests/test_estimate_a.py index 5ad6f43..8a23b4b 100644 --- a/seismostats/analysis/tests/test_estimate_a.py +++ b/seismostats/analysis/tests/test_estimate_a.py @@ -1,6 +1,7 @@ import numpy as np import warnings from numpy.testing import assert_almost_equal +from datetime import datetime, timedelta from seismostats.analysis.estimate_a import ( estimate_a_classic, estimate_a_positive) @@ -32,14 +33,19 @@ def test_estimate_a_classic(): def test_estimate_a_positive(): - mags = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) - times = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) + mags = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 10]) + # array of dates from 1.1.2000 to 10.1.2000 + times = np.arange(datetime(2000, 1, 1), datetime( + 2000, 1, 12), timedelta(days=1)).astype(datetime) + + print(times) + print(type(times)) a = estimate_a_positive(mags, times, delta_m=1) - assert_almost_equal(10**a, 9.0) + assert_almost_equal(10**a, 10.0) # reference magnitude is given and b-value given a = estimate_a_positive(mags, times, delta_m=1, mc=1, m_ref=0, b_value=1) - assert_almost_equal(10**a, 90.0) + assert_almost_equal(10**a, 100.0) # reference magnitude but no b-value try: @@ -49,7 +55,7 @@ def test_estimate_a_positive(): # reference time is given a = estimate_a_positive(mags, times, delta_m=1, scaling_factor=10) - assert_almost_equal(10**a, 0.9) + assert_almost_equal(10**a, 1.0) # magnitudes not cut at mc with warnings.catch_warnings(record=True) as w: From 2f5991a3a69cbd016978df59fcaf3d51ac07d9dc Mon Sep 17 00:00:00 2001 From: Aron Mirwald Date: Mon, 19 Aug 2024 20:07:14 +0200 Subject: [PATCH 8/8] delete print --- seismostats/analysis/tests/test_estimate_a.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/seismostats/analysis/tests/test_estimate_a.py b/seismostats/analysis/tests/test_estimate_a.py index 8a23b4b..7fdd222 100644 --- a/seismostats/analysis/tests/test_estimate_a.py +++ b/seismostats/analysis/tests/test_estimate_a.py @@ -38,8 +38,6 @@ def test_estimate_a_positive(): times = np.arange(datetime(2000, 1, 1), datetime( 2000, 1, 12), timedelta(days=1)).astype(datetime) - print(times) - print(type(times)) a = estimate_a_positive(mags, times, delta_m=1) assert_almost_equal(10**a, 10.0)