On this page I'm providing Python scripts to simulate and fit strain and seismicity data associated with "acceleration-to-failure" processes. If you plan to use these scripts as part of publication, please reference Bell et al. 2011 (GRL).

* NEW!!!* - Updated in June 2013 to include maximum-likelihood forecasting methods from Bell et al. 2013 (GJI). here

Here are two Python scripts to generate:

- Earthquake or AE data as a heterogeneous Poisson process, with a power-law or exponential acceleration in mean rate and with magnitudes distributed according to the Gutenberg-Richter distribution
- Strain data with a Gaussian measurement uncertainty, with a power-law acceleration in mean rate

Two example data files are also provided.

- Python script for generating and plotting synthetic earthquake/AE data
- Python script for generating and plotting synthetic strain data

Here are Python scripts to fit power-law accelerations to failure in strain and earthquake/AE data using different methods (parameter confidence intervals to follow):

- To strain data using a GLM routine with log-link function and Gaussian errors (uses scikits.statsmodels)
- To strain data by minimizing the least-squares resdiuals between model and data (uses scipy.optimize)
- To earthquake/AE data using a GLM routine with log-link function and Poisson errors (uses scikits.statsmodels)
- To earthquake/AE data using by minimizing Ogata's log-likelihood function (uses scipy.optimize)

- Python script for fitting power-law to strain using GLM
- Python script for fitting power-law to strain by minimizing LS
- Python script for fitting power-law to earthquake rates using GLM
- Python script for fitting power-law to earthquake occurrence times using Ogata's method

Information criterion are a pragmatic tool for comparing the performance of models, taking model complexity into account. Here I provide a code for using the Akaike information criterion (AIC) for comparing a power-law and exponential model for accelerating rates of earthquakes (using Ogata's maximum likelihood method to fit the models). The preferred model is the one with the lower AIC; see Bell et al. 2011 (GJI) for examples. Note that rigorous model comparison requires calculation of the Bayesian "evidence" (see Toussaint, 2011, for discussion).

Bell et al. 2013 (GJI) outlines maximum-likelihood methods for fitting and forecasting accelerating rates of earthquakes.

These scripts are written in Python. Python is a general-purpose, readable, high-level programming language. It is open source, with community development. NumPy is an extension to Python that allows operation on multi-dimensional arrays and matrices. SciPy is a library of mathematical tools. Matplotlib is a plotting library. Together these components provide a free, open source alternative to Matlab.

For windows users I recommend the Pythonxy distribution of Python, based on Qt and Spyder.

Some of these scripts use the Statsmodels module, and the Optimize module of SciPy.

- Bell, A.F., M. Naylor, M.J. Heap and I.G. Main,
*Forecasting volcanic eruptions and other material failure phenomena: an evaluation of the failure forecast method*, Geophys. Res. Lett., 38, L15304, doi:10.1029/2011GL048155**(2011)**. - Bell, A.F., J. Greenhough, M.J. Heap and I.G. Main,
*Challenges for forecasting based on accelerating rates of earthquakes at volcanoes and laboratory analogues*, Geophys. J. Int., 185, 2, doi:10.1111/j.1365-246X.2011.04982.x,**(2011)**. - Toussaint, U. v.,
*Bayesian inference in physics*, Rev. Mod. Phys., 83, doi:10.1103/RevModPhys.83.943,**(2011)**. - Ogata, Y.,
*Statistical analysis of seismicity: updated version (SASeis2006)*, Computer Science Monographs, 33**(2006)**.