swiss-seismological-service / SeismoStats

https://seismostats.readthedocs.io
GNU Affero General Public License v3.0
7 stars 3 forks source link

Possible(?) minor bug in b-positive #158

Closed thpap closed 4 months ago

thpap commented 4 months ago

Hi guys,

Line 289 of estimate_beta.py: mag_diffs = abs(mag_diffs[mag_diffs > 0])

Shouldn't the inequality for binned magnitudes be over the binning width and not over 0? See text below equation (8) and Appendix of Elst (2021).

I suggest changing this to: mag_diffs = abs(mag_diffs[(mag_diffs > delta_m) | (np.isclose(mag_diffs, delta_m, atol = 0.001))])

The "or" operator because floating numbers mess it up if you use ">=".

aronsho commented 4 months ago

I agree that it should be the binning width- in most cases, this is also what the code does: if the magnitudes are binned correctly, mag_diffs >0 results in the same as mag_diffs >= delta_m.

Maybe, we could put mag_diffs > delta_m/2, this would make sure that some mean rounding errors also would be excluded..

thpap commented 4 months ago

Hmm ok yes I see what you mean. Maybe adding a small tolerance like you suggested, or catching it at the upper bound with an "or" logical operator is more mistake-proof, but you are right :) Thanks for checking!