linksrm: Linked Stress Release Model Object
In PtProcess: Time Dependent Point Process Modelling

Description Usage Arguments Details Examples

Creates a point process model object with class "linksrm".

1	linksrm(data, gif, marks, params, gmap, mmap, TT)

`data`	a `data.frame` containing the history of the process, denoted below as Ht. It should contain all variables that are required to evaluate the `gif` function and the mark distribution, though can contain others too. No history is represented as `NULL`.
`gif`	ground intensity function. At this stage, this can only be `linksrm_gif` or modifications of that function; see “Details” below.
`marks`	mark distribution. See topic `marks` for further details.
`params`	numeric vector of all model parameters.
`gmap`	`expression`, maps the model parameters (`params`) into the parameter sub-space of the ground intensity function; see “Details” below.
`mmap`	`expression`, maps the model parameters (`params`) into the parameter sub-space of the mark distribution; see “Details” below.
`TT`	vector of length 2, being the time interval over which the integral of the ground intensity function is to be evaluated.

The linked stress release model has a slightly peculiar structure which makes it difficult to fit into the mpp class. While the region should be thought of as a mark, it is completely defined by the function linksrm_gif, and hence from the programming perspective the region mark is really tied in with the gif function. Hence at the moment, the linked stress release model is treated as a special case. There may be other models that could be grouped into this class.

p <- c(-1.5, -1.5, 0.01, 0.03, 2, -0.5, 0.2, 1, 1*log(10), 3)
TT <- c(0, 1000)

rexptrunc_mark <- function(ti, data, params){
    x <- rexp(n=1, params[1])
    x[x > params[2]] <- params[2]
    names(x) <- "magnitude"
    return(x)
}

x <- linksrm(data=NULL,
             gif=linksrm_gif,
             marks=list(NULL, rexptrunc_mark),
             params=p,
             gmap=expression(params[1:8]),
             mmap=expression(params[9:10]),
             TT=TT)

x <- simulate(x, seed=5)
print(logLik(x))

#   estimate parameters
temp_map <- function(y, p){
    #    map only gif parameters into model object
    y$params[1:8] <- p
    return(y)
}

weight <- c(0.1, 0.1, 0.005, 0.005, 0.1, 0.1, 0.1, 0.1)

#   see manual page for linksrm_gif for modifications to
#   make calculations faster

#   for testing, restrict to 5 iterations
z <- nlm(neglogLik, p[1:8], object=x, pmap=temp_map,
         hessian=TRUE, gradtol=1e-08, steptol=1e-10,
         print.level=2, iterlim=5, typsize=weight)

param.names <- c("a1", "a2", "b1", "b2", "c11", "c12", "c21", "c22")
param.est <- cbind(p[1:8], z$estimate, sqrt(diag(solve(z$hessian))))
dimnames(param.est) <- list(param.names,
                            c("Actual", "Estimate", "StdErr"))
print(param.est)

#   place parameter estimates into model object
x <- temp_map(x, z$estimate)

#   plot ground intensity function
par.default <- par(mfrow=c(2,1), mar=c(4.1, 4.1, 0.5, 1))
x$gif <- linksrm_gif
plot(x, 1, xlab="")
plot(x, 2)
par(par.default)

#   plot "residuals" for each region
tau <- residuals(x)
par(mfrow=c(2,1))
for (i in 1:2){
    plot(tau[[i]], ylab="Transformed Time",
         xlab="Event Number", main=paste("Region", i))
    abline(a=0, b=1, lty=2, col="red")
}

#   plot cusum of "residuals" for each region
for (i in 1:2){
    plot(tau[[i]] - 1:length(tau[[i]]), ylab="Cusum of Transformed Time",
         xlab="Event Number", main=paste("Region", i))
    abline(h=0, lty=2, col="red")
}

par(mfrow=c(1,1))

[1] -1136.044
iteration = 0
Step:
[1] 0 0 0 0 0 0 0 0
Parameter:
[1] -1.50 -1.50  0.01  0.03  2.00 -0.50  0.20  1.00
Function Value
[1] 1136.044
Gradient:
[1]    3.140403   -5.617850 -580.972142   43.155440  -36.387405  -38.165589
[7]  180.092309  222.638991

iteration = 1
Step:
[1] -3.448537e-05  6.169070e-05  1.594942e-05 -1.184746e-06  3.995771e-04
[6]  4.191037e-04 -1.977629e-03 -2.444842e-03
Parameter:
[1] -1.50003449 -1.49993831  0.01001595  0.02999882  2.00039958 -0.49958090
[7]  0.19802237  0.99755516
Function Value
[1] 1135.576
Gradient:
[1]    2.406780    7.629559 -587.210124  134.937106  -32.953252  -33.999518
[7]  -14.957625   -9.464118

iteration = 2
Step:
[1] -2.687337e-05 -7.859705e-05  1.619549e-05 -3.626006e-06  3.651475e-04
[6]  3.770046e-04  8.251514e-05  6.200100e-06
Parameter:
[1] -1.50006136 -1.50001691  0.01003214  0.02999519  2.00076472 -0.49920389
[7]  0.19810489  0.99756136
Function Value
[1] 1135.54
Gradient:
[1]    1.746726    7.345193 -592.048043  131.962111  -29.855719  -30.242831
[7]  -10.907108   -4.654837

iteration = 3
Step:
[1] -0.0002993517 -0.0011088792  0.0002370780 -0.0000517776  0.0048792376
[6]  0.0049604497  0.0007242595 -0.0004737770
Parameter:
[1] -1.50036071 -1.50112579  0.01026922  0.02994341  2.00564396 -0.49424344
[7]  0.19882915  0.99708758
Function Value
[1] 1135.308
Gradient:
[1]   -6.9285415    6.4416582 -605.7352783  118.9270912   11.7361724
[6]   20.1975861    0.2160682    8.4291700

iteration = 4
Step:
[1]  6.018383e-05 -1.131038e-04  2.527649e-05 -5.212450e-06  7.363099e-05
[6] -1.138430e-05  4.664415e-05 -8.342591e-05
Parameter:
[1] -1.5003005 -1.5012389  0.0102945  0.0299382  2.0057176 -0.4942548  0.1988758
[8]  0.9970042
Function Value
[1] 1135.291
Gradient:
[1]   -6.944554    6.542231 -602.156605  119.350267   11.918747   20.403514
[7]   -1.495183    6.377339

iteration = 5
Parameter:
[1] -1.49166415 -1.51660782  0.01370565  0.02923372  2.01318178 -0.49880169
[7]  0.20662805  0.98746575
Function Value
[1] 1133.909
Gradient:
[1]   -3.725711   10.221956 -114.047181   97.102176   10.771722   15.335716
[7]  -84.403417  -94.348895

Iteration limit exceeded.  Algorithm failed.

    Actual    Estimate      StdErr
a1   -1.50 -1.49166415 0.139289690
a2   -1.50 -1.51660782 0.121134114
b1    0.01  0.01370565 0.004794368
b2    0.03  0.02923372 0.005478861
c11   2.00  2.01318178 0.411117632
c12  -0.50 -0.49880169 0.334259953
c21   0.20  0.20662805 0.190716826
c22   1.00  0.98746575 0.158255452