Calculates the value of the ground intensity of a Linked Stress Release Model (LSRM). This model allows for multiple linked regions, where the stress can be transferred between the regions.
linksrm_gif(data, evalpts, params, TT=NA, tplus=FALSE, eta=0.75)a data frame containing the event history, where each row represents one event. There must be columns named "time", usually the number of days from some origin; "magnitude" which is the event magnitude less the magnitude threshold, i.e. \(M_k - M_0\); and "region" which are consecutively numbered starting at 1.
a matrix or data.frame. It must include two columns named "time" and "region" that can be referred to as evalpts[,"time"] and evalpts[,"region"], respectively. The function will be evaluated at these points.
vector of parameters of length \(n^2+2n\), where \(n\) is the number of regions, for the proposed LSRM in the following order: $$ (a_1, \cdots, a_n, b_1, \cdots, b_n, c_{11}, c_{12}, c_{13}, \cdots, c_{nn}). $$
vector of length 2, being the time interval over which the integral of the ground intensity function is to be evaluated.
logical, \(\lambda_g(t,i|{\cal H}_t)\) is evaluated as \(\lambda_g(t^+,i|{\cal H}_t)\) if TRUE, else \(\lambda_g(t^-,i|{\cal H}_t)\).
a scalar used in the stress calculations, see Details below.
Two usages are as follows.
linksrm_gif(data, evalpts, params, tplus=FALSE, eta=0.75) linksrm_gif(data, evalpts=NULL, params, TT, eta=0.75)
The first usage returns a vector containing the values of \(\lambda_g(t,i)\) evaluated at the specified “time-region” points. In the second usage, it returns a vector containing the value of the integral for each region.
rateis "increasing".
regionsis expression(sqrt(length(params) + 1) - 1).
The function linksrm_gif calculates the stress reduction matrices St1 and St2 every time that the function is called. Ideally, these should be calculated once and be included within the model object. Currently, the structure of the model object is not sufficiently flexible. However, the user can create a new function to calculate St1 and St2 once. This will only work if the event history is not changing between successive calls (e.g. parameter estimation). However, in a simulation, the history changes with the addition of each new event, and in this situation St1 and St2 need to be calculated with every function call.
The modified function, as described below, will write the objects St1 and St2 to a temporary database (position 2 in the search path). Consequently, it cannot be defined within the package itself because this violates the CRAN rules. The function linksrm_gif contains markers indicating the beginning and ending of the parts where St1 and St2 are calculated. The modified function is made by editing the function linksrm_gif. We firstly deparse the function linksrm_gif (i.e. put the contents into a character vector). We initially create a temporary database called PtProcess.tmp in which to write St1 and St2. We then search for the line numbers that mark the beginning and ending of the parts where St1 and St2 are calculated. We replace the beginning of each with a conditional statement so that the contents are only run if these two objects do not already exist. We then parse the lines of code in the character vector back into a function, and call this new function linksrm1_gif. The same thing can be achieved by dumping linksrm_gif to a text file and editing manually.
# define linksrm1_gif by modifying linksrm_gif# put function linksrm_gif into a character vector tmp <- deparse(linksrm_gif)
# remove "if (FALSE)" lines linenum <- grep("if \(FALSE\)", tmp) tmp <- tmp[-linenum]
# attach new database at pos=2 in search path called PtProcess.tmp linenum <- grep("attach new database to search path", tmp) tmp[linenum] <- "if (!any(search()==\"PtProcess.tmp\")) attach(NULL, pos=2L, name=\"PtProcess.tmp\", warn.conflicts=TRUE)"
# calc St1 if St1 does not exist linenum <- grep("this loop calculates St1", tmp) tmp[linenum] <- "if (!exists(\"St1\", mode = \"numeric\")) {" linenum <- grep("assign statement for St1", tmp) tmp[linenum] <- "assign(\"St1\", St1, pos=\"PtProcess.tmp\")" linenum <- grep("end loop St1", tmp) tmp[linenum] <- "}"
# calc St2 if St2 does not exist linenum <- grep("this loop calculates St2", tmp) tmp[linenum] <- "if (!exists(\"St2\", mode = \"numeric\")) {" linenum <- grep("assign statement for St2", tmp) tmp[linenum] <- "assign(\"St2\", St2, pos=\"PtProcess.tmp\")" linenum <- grep("end loop St2", tmp) tmp[linenum] <- "}"
linksrm1_gif <- eval(parse(text=tmp))
Warning: The function linksrm1_gif checks to see whether the matrices St1 and St2 exist. If so, these existing matrices are used, and new ones are not calculated. Therefore when using linksrm1_gif for parameter estimation, one must check for the existence of such matrices, and delete upon starting to fit a new model:
if (exists("St1")) rm(St1)
if (exists("St2")) rm(St2)
or detach the database as detach(2). The objects St1 and St2 will exist for the duration of the current R session, so should be deleted when no longer required.
The ground intensity for the \(i\)th region is assumed to have the form
$$
\lambda_g(t,i | {\cal H}_t) = \exp\left\{ a_i + b_i\left[t - \sum_{j=1}^n c_{ij} S_j(t)\right]\right\}
$$
with params\(= \code{c}(a_1, \cdots, a_n, b_1, \cdots, b_n, c_{11}, c_{12}, c_{13}, \cdots, c_{nn})\); and
$$
S_j(t) = \sum_k 10^{\eta(M_k-M_0)},
$$
where the summation is taken over those events in region \(j\) with time \(t_k < t\). This model has been discussed by Bebbington & Harte (2001, 2003). The default value of \(\eta = \code{eta} = 0.75\).
Cited references are listed on the PtProcess manual page.
General details about the structure of ground intensity functions are given in the topic gif.