Description Usage Arguments Details Examples
Creates a point process model object with class "linksrm"
.
1 |
data |
a |
gif |
ground intensity function. At this stage, this can only be |
marks |
mark distribution. See topic |
params |
numeric vector of all model parameters. |
gmap |
|
mmap |
|
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 | 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.