Source code for openquake.hmtk.seismicity.occurrence.aki_maximum_likelihood
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## LICENSE## Copyright (C) 2010-2023 GEM Foundation, G. Weatherill, M. Pagani,# D. Monelli.## The Hazard Modeller's Toolkit is free software: you can redistribute# it and/or modify it under the terms of the GNU Affero General Public# License as published by the Free Software Foundation, either version# 3 of the License, or (at your option) any later version.## You should have received a copy of the GNU Affero General Public License# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>## DISCLAIMER## The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein# is released as a prototype implementation on behalf of# scientists and engineers working within the GEM Foundation (Global# Earthquake Model).## It is distributed for the purpose of open collaboration and in the# hope that it will be useful to the scientific, engineering, disaster# risk and software design communities.## The software is NOT distributed as part of GEM’s OpenQuake suite# (https://www.globalquakemodel.org/tools-products) and must be considered as a# separate entity. The software provided herein is designed and implemented# by scientific staff. It is not developed to the design standards, nor# subject to same level of critical review by professional software# developers, as GEM’s OpenQuake software suite.## Feedback and contribution to the software is welcome, and can be# directed to the hazard scientific staff of the GEM Model Facility# (hazard@globalquakemodel.org).## The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed WITHOUT# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License# for more details.## The GEM Foundation, and the authors of the software, assume no# liability for use of the software.importwarningsimportnumpyasnpfromopenquake.hmtk.seismicity.occurrence.baseimport(SeismicityOccurrence,OCCURRENCE_METHODS,)fromopenquake.hmtk.seismicity.occurrence.utilsimport(recurrence_table,input_checks,)
[docs]defcalculate(self,catalogue,config=None,completeness=None):""" Calculation of b-value and its uncertainty for a given catalogue, using the maximum likelihood method of Aki (1965), with a correction for discrete bin width (Bender, 1983). :param catalogue: See :class:`openquake.hmtk.seismicity.occurrence.base.py` for further explanation :param config: The configuration in this case do not contains specific information :keyword float completeness: Completeness magnitude :return float bval: b-value of the Gutenberg-Richter relationship :return float sigma_b: Standard deviation of the GR b-value """# Input checks_cmag,_ctime,_ref_mag,dmag,config=input_checks(catalogue,config,completeness)rt=recurrence_table(catalogue.data["magnitude"],dmag,catalogue.data["year"])bval,sigma_b=self._aki_ml(rt[:,0],rt[:,1])returnbval,sigma_b
def_aki_ml(self,mval,number_obs,dmag=0.1,m_c=0.0):""" :param numpy.ndarray mval: array of reference magnitudes (column 0 from recurrence table) :param numpy.ndarray number_obs: number of observations in magnitude bin (column 1 from recurrence table) :keyword float dmag: magnitude interval :keyword float m_c: completeness magnitude :return float bval: b-value of the Gutenberg-Richter relationship :return float sigma_b: Standard deviation of the GR b-value """# Exclude data below Mcid0=mval>=m_cmval=mval[id0]number_obs=number_obs[id0]# Get Number of events, minimum magnitude and mean magnitudeneq=np.sum(number_obs)ifneq<=1:# Cannot determine b-value (too few event) return NaNswarnings.warn("Too few events (<= 1) to calculate b-value")returnnp.nan,np.nanm_min=np.min(mval)m_ave=np.sum(mval*number_obs)/neq# Calculate b-valuebval=np.log10(np.exp(1.0))/(m_ave-m_min+(dmag/2.0))# Calculate sigma b from Bender estimatorsigma_b=np.sum(number_obs*((mval-m_ave)**2.0))/(neq*(neq-1))sigma_b=np.log(10.0)*(bval**2.0)*np.sqrt(sigma_b)returnbval,sigma_b