simulate: Simulate a Point Process In PtProcess: Time Dependent Point Process Modelling

Description

Provides methods for the generic function `simulate`.

Usage

 ```1 2 3 4 5 6``` ```## S3 method for class 'mpp' simulate(object, nsim=1, seed=NULL, max.rate=NA, stop.condition=NULL, ...) ## S3 method for class 'linksrm' simulate(object, nsim=1, seed=NULL, max.rate=NA, stop.condition=NULL, ...) ```

Arguments

 `object` an object with class `"mpp"` or `"linksrm"`. `nsim` has no effect, and is only included for compatibility with the generic function `simulate`. See section “Length of Simulated Series” below for control information. `seed` seed for the random number generator. `max.rate` maximum rate, only used if the attribute of `object\$gif` is `"bounded"`. It is the maximum value of `object\$gif` on the simulation interval `object\$TT`. `stop.condition` a function returning a logical value. It is called after the addition of each simulated event. The simulation continues until either `object\$TT[2]` is exceeded or `stopping.condition` returns `TRUE`. See section “Length of Simulated Series” below for further information. `...` other arguments.

Details

The thinning method (Ogata, 1981; Lewis & Shedler, 1979) is used to simulate a point process with specified ground intensity function. The method involves calculating an upper bound for the intensity function, simulating a value for the time to the next possible event using a rate equal to this upper bound, and then calculating the intensity at this simulated point; hence these “events” are simulated too frequently. The ratio of this rate with the upper bound is compared with a uniform random number to randomly determine whether the simulated time is retained or not (i.e. thinned).

The functions need to calculate an upper bound for the intensity function. The ground intensity functions will usually be discontinuous at event times, but may be monotonically increasing or decreasing at other times. The ground intensity functions have an attribute called `rate` with values of `"bounded"`, `"increasing"` or `"decreasing"`. This information is used to determine the required upper bounded.

The function `simulate.linksrm` is currently only used in conjunction with `linksrm_gif`, or a variation of that function. It expects the `gif` function to have an attribute called `regions`, which may be an expression, being the number of regions. The method used by the function `simulate.linksrm` also assumes that the function is “increasing” (i.e. rate, summed over all regions, apart from discontinuous jumps), hence a positive tectonic input over the whole system.

Value

The returned value is an object of the same class as `object`. It will contain all events prior to `object\$TT[1]` in `object\$data` and all subsequently simulated events. Variables (columns) in `object\$data` will be restricted to `"time"` and those for which a mark is simulated.

Length of Simulated Series

The interval of time over which events are simulated is determined by `object\$TT`. Simulation starts at `object\$TT[1]` and stops at `object\$TT[2]`. The “current” dataset will consist of all events prior to `object\$TT[1]` in `object`, plus subsequently simulated events. A more complicated stopping condition can be formulated by using the argument `stop.condition`.

The argument `stop.condition` can be assigned a function that returns a logical value. The assigned function is a function of the “current” dataset. It is executed near the bottom of `simulate.mpp` (check by printing the function). Simulation will then continue until either the stopping condition has been met or the current time exceeds `object\$TT[2]`.

For example, we may want to simulate until the first earthquake with a magnitude of 8. Assume that the current dataset contains a variable with name `"magnitude"` (untransformed). We would then assign `Inf` to `object\$TT[2]`, and write this condition as a function:

 ```1 2 3 4 5``` ``` stop.cond <- function(data){ n <- nrow(data) # most recent event is the nth return(data\$magnitude[n] >= 8) } ```

References

Cited references are listed on the PtProcess manual page.

Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20``` ```TT <- c(0, 1000) bvalue <- 1 params <- c(-2.5, 0.01, 0.8, bvalue*log(10)) x <- mpp(data=NULL, gif=srm_gif, marks=list(NULL, rexp_mark), params=params, gmap=expression(params[1:3]), mmap=expression(params[4]), TT=TT) x <- simulate(x, seed=5) y <- hist(x\$data\$magnitude, xlab="Magnitude", main="") # overlay with an exponential density magn <- seq(0, 3, length.out=100) points(magn, nrow(x\$data)*(y\$breaks[2]-y\$breaks[1])* dexp(magn, rate=1/mean(x\$data\$magnitude)), col="red", type="l") ```

Example output

```
```

PtProcess documentation built on Nov. 17, 2017, 7:12 a.m.