Help with the Weichert occurrence algorithm

[Edwar Q S]

Help with the Weichert occurrence algorithm
Good morning. I am trying to use the Weichert occurrence algorithm but I have noticed that the value of 'a' is extremely high. Is this due to the completeness table or what could it be due to? I attach the script I am trying to use, a capture of the error associated with the high value of 'a' and the csv file with my catalog. If the reason is the completeness table, could you guide me how to define this in the best way please? Thank you so much.

Python code:

from openquake.hmtk.seismicity.catalogue import Catalogue
from openquake.hmtk.seismicity.completeness.comp_stepp_1971 import Stepp1971
from openquake.hmtk.plotting.seismicity.completeness.plot_stepp_1972 import create_stepp_plot

import numpy as np
import os
from openquake.hmtk.parsers.catalogue.csv_catalogue_parser import CsvCatalogueParser
import matplotlib.pyplot as plt
from cycler import cycler

# Seismicity tools: Recurrence methods
from openquake.hazardlib.mfd import TruncatedGRMFD
from openquake.hmtk.seismicity.occurrence.aki_maximum_likelihood import AkiMaxLikelihood
from openquake.hmtk.seismicity.occurrence.b_maximum_likelihood import BMaxLikelihood
from openquake.hmtk.seismicity.occurrence.kijko_smit import KijkoSmit
from openquake.hmtk.seismicity.max_magnitude.kijko_sellevol_bayes import KijkoSellevolBayes
from openquake.hmtk.seismicity.occurrence.weichert import Weichert
from openquake.hmtk.plotting.seismicity.occurrence.recurrence_plot import plot_recurrence_model

# Path to the CSV file
csv_file_path = 'Filtrado_FS2_4.5.csv'

# Create an instance of the CSV catalogue parser
parser = CsvCatalogueParser(input_file=csv_file_path)

# Read the CSV file and obtain a Catalogue object
catalogue = parser.read_file()

# Check the number of events in the catalogue
print("The catalogue contains %g events" % catalogue.get_number_events())

stepp = Stepp1971()

# Configure the completeness analysis parameters
completeness_config = {
    'magnitude_bin': 0.5,
    'time_bin': 5,
    'increment_lock': True
}

# Execute the completeness analysis using Stepp (1971)
print('Executing Stepp (1971) completeness analysis:')
completeness_table = stepp.completeness(catalogue, completeness_config)
print(completeness_table)
print('Done!')

# completeness_table=np.array([[1982. , 4.5] ,[1961. , 5.0], [1945. , 5.5] , [1919. , 6.], [1912. , 6.5], [1910. , 7.0] ])

recurrence_estimator = Weichert()
recurrence_config = {"magnitude_interval": 0.1}

bval, sigma_b, aval, sigma_a = recurrence_estimator.calculate(catalogue,
                                                              recurrence_config,
                                                              completeness_table)

print("a = %.3f (+/- %.3f),  b = %.3f (+/-%.3f)" % (aval, sigma_a, bval, sigma_b))

mfd0 = TruncatedGRMFD(4.5, 7.6, 0.1, aval, bval)

[Peter Pažák]

Hi, I do not think the reason is the completeness table, BMaxLikelihood method with the same completeness_table gathers  a = 7.004 (+/- 0.125), b = 1.164 (+/-0.027),

which seems to fit the data quite reasonably. My guess is the Wiechert algorithm is not converging for the data, not sure how this could be improved - maybe it is better to just use those other algorithms?

Peter

Dátum: štvrtok 13. júna 2024, čas: 10:05:08 UTC+2, odosielateľ: edwa...@gmail.com

[kendra....@globalquakemodel.org]

A while ago there was a change to how this algorithm works. 

The part recurrence_estimator.calculate must be changed to recurrence_estimator.calc

Using this with the rest of your code, I get a=8.506