Implement method for Mc by b-value stability (based on Cao & Gao 2002). #165
Conversation
seismostats/analysis/estimate_mc.py
Outdated
| Args: | ||
| sample : np.array Vector of magnitudes | ||
| delta_m : float. Discretization of the magnitudes. | ||
| stability_factor : float. Magnitude bin to consider for the |
There was a problem hiding this comment.
I would prefer to use a stability_range- and the binning would be used as is given by the delta_m.
The stability_range would be default at 0.5, but could be any value, (e.g., 0.2), and stability_range/delta_m points will be evaluated for stability. E.g., for stability_range = 0.6 and delta_m = 0.1, six points will be evaluated, but if delta_m= 0.2, only 3 points.
This would make the parameter stability_range more meaningful. What do you think?
There was a problem hiding this comment.
Would also be very easy to implement I think :)
There was a problem hiding this comment.
Yes I agree it would be more elegant and more in line with what we decided last time (i.e. not sticking to the original Cao & Gao 2002 definition).
I'll modify it accordingly.
| for each step in Mc | ||
| bs : np.ndarray. Estimated b-value for each step in Mc | ||
| """ | ||
| sample = bin_to_precision(sample, delta_x=delta_m) |
There was a problem hiding this comment.
I think it would be better to test if it is binned correctly and then give a warning. Maybe this is for the future though, because this would require a test_binning function
There was a problem hiding this comment.
AFAIK this is inconsistent between methods. Some require the catalogue to be pre-binned-to-precision and some don't. I would keep this out for a different clean-up issue.
seismostats/analysis/estimate_mc.py
Outdated
| b = estimate_b(sample[sample >= mcp - delta_m / 2], | ||
| mcp, delta_m, b_parameter='b_value', method="tinti") | ||
| b_ex.append(b) | ||
| b_avg = np.sum(b_ex) / 5 |
There was a problem hiding this comment.
Here, it should be divided by len(b_ex)? It seems that it is possible that there are a different number of b-values than 5
There was a problem hiding this comment.
true, in line with our new and updated definition that doesn't stick to the original definition.
|
|
||
|
|
||
| @pytest.fixture | ||
| def swiss_2023_magnitudes(): |
There was a problem hiding this comment.
I don't like that we are dependent that the server is always working for the testing. I We should instead have a static catalog (I guess a csv?)
There was a problem hiding this comment.
To be discussed with Nicolas at our next meeting.
| @@ -207,3 +236,11 @@ def test_estimate_mc_maxc(setup_magnitudes): | |||
| mc = mc_max_curvature(setup_magnitudes, delta_m=0.1, correction_factor=0.2) | |||
|
|
|||
| assert_equal(1.3, mc) | |||
|
|
|||
|
|
|||
| def test_estimate_mc_bvalue_stability(swiss_2023_magnitudes): | |||
There was a problem hiding this comment.
This should be tested also for different binning and stability range :)
There was a problem hiding this comment.
Another missing test: one catalogue that fails to be stable enough.
There was a problem hiding this comment.
test for 0.1 binning added
…logical-service/SeismoStats into mc_bvalue_stability
|
@aronsho should be good to go now (minus the test for failing catalogue but we can figure this one out later) |
The stability test is done on half a magnitude unit and not 5*magnitude binning.
Also new: introduced pytest fixture for future use of the swiss 2023 catalogue in tests for Mc methods (GOF I'm coming for you...)