library("PtProcess") data("Phuket") # magnitudes are rounded to 1 dp # threshold magnitude = M0 = 4.95 Phuket$magnitude <- Phuket$magnitude - 4.95 #---------------------------------------------------------- # Define Full Model Object dmagn_mark <- function(x, data, params){ # Gamma distribution # exponential density when params[7]=0 if (params[7]>0){ lambda <- etas_gif(data, x[,"time"], params=params[1:5]) y <- dgamma(x[,"magnitude"], shape=1+sqrt(lambda)*params[7], rate=params[6], log=TRUE) } else y <- dexp(x[,"magnitude"], rate=params[6], log=TRUE) return(y) } rmagn_mark <- function(ti, data, params){ # Gamma distribution # exponential density when params[7]=0 if (params[7]>0){ lambda <- etas_gif(data, ti, params=params[1:5]) y <- rgamma(1, shape=1+sqrt(lambda)*params[7], rate=params[6]) } else y <- rexp(1, rate=params[6]) return(list(magnitude=y)) } TT <- c(0, 1827) params <- c(0.05, 5.2, 1.8, 0.02, 1.1, 1/mean(Phuket$magnitude), 0) x <- mpp(data=Phuket, gif=etas_gif, mark=list(dmagn_mark, rmagn_mark), params=params, TT=TT, gmap=expression(params[1:5]), mmap=expression(params)) #---------------------------------------------------------- # Fit Model with Exponential Magnitude Distribution expmap <- function(y, p){ # for exponential distribution y$params[1:5] <- exp(p) return(y) } initial <- log(params[1:5]) z <- optim(initial, neglogLik, object=x, pmap=expmap, control=list(trace=1, maxit=100)) initial <- z$par z <- nlm(neglogLik, initial, object=x, pmap=expmap, print.level=2, iterlim=500, typsize=initial) # write estimates to a new model object x0 x0 <- expmap(x, z$estimate) print(logLik(x0)) #---------------------------------------------------------- # Fit Model with Gamma Magnitude Distribution allmap <- function(y, p){ y$params <- exp(p) return(y) } initial <- log(c(0.05, 5.2, 1.8, 0.02, 1.1, 1/mean(Phuket$magnitude), 0.1)) z <- optim(initial, neglogLik, object=x, pmap=allmap, control=list(trace=1, maxit=200)) initial <- z$par z <- nlm(neglogLik, initial, object=x, pmap=allmap, print.level=2, iterlim=500, typsize=initial) # write estimates to a new model object x1 x1 <- allmap(x, z$estimate) print(logLik(x1)) #---------------------------------------------------------- # Print Summary of Both Models print(summary(x0)) print(summary(x1)) #---------------------------------------------------------- # Table of Parameter Estimates param.est <- rbind(cbind(x0$params, x1$params), c(logLik(x0), logLik(x1))) colnames(param.est) <- c("Exp Model", "Gamma Model") rownames(param.est) <- c("p1=mu", "p2=A", "p3=alpha", "p4=c", "p5=p", "p6", "p7", "logLik") print(param.est) #---------------------------------------------------------- # Simulations # Use fitted model in oject x2 (gamma case) x2 <- x1 # Define simulation interval # from 1 Jan 2009 (x1$TT[2]) to infinity x2$TT <- c(x2$TT[2], Inf) stop.cond <- function(data){ # stop if the last simulated event has # a magnitude >= 6.5 (4.95+1.55) n <- nrow(data) return(data$magnitude[n] >= 1.55) } # y[i] is the time to the next event in the ith simulation # that satisfies the stopping condition y <- rep(NA, 2000) for (i in 1:2000){ print(i) x3 <- simulate(x2, seed=i, stop.cond=stop.cond) y[i] <- x3$data$time[nrow(x3$data)] } #---------------------------------------------------------- # Generate Graphics library("mapdata") #---------------------------------------------------------- # Define Big Events big <- list() big$id <- c(14, 35, 511, 733, 875, 972, 1124, 1194) big$time <- c(206.6079, 360.0409, 452.6733, 570.6542, 866.6447, 1350.4656, 1516.3587, 1639.4863) big$date <- c("25Jul04", "26Dec04", "28Mar05", "24Jul05", "16May06", "12Sep07", "20Feb08", "27Jun08") #---------------------------------------------------------- # Magnitude-time & GIF plots pdf("phuket1.pdf", height=9.5, width=8) par(mfrow=c(2,1), mar=c(4.1, 4.1, 0, 0.2), oma=c(0, 0, 4.5, 0)) plot(x0$data$time, x0$data$magnitude+4.95, type="h", xlim=x0$TT, xlab="", ylab="Event Magnitude") axis(3, at=big$time, labels=big$date, cex.axis=1, las=2) plot(x0, log=TRUE, xlim=x0$TT, xlab="t (Days Since 1 Jan 2004)") dev.off() #---------------------------------------------------------- # Residual Process Std Plot pdf(file="phuket2.pdf", height=9, width=8) # note that residuals(x1) returns a "ts" object par(mar=c(4, 4, 5, 0.2)) plot(residuals(x1), xlab="Event Number", ylab="Transformed Time", pty="s") points(residuals(x0), lty=2, type="l", col="red") axis(3, at=big$id[-1], labels=big$date[-1], las=2) abline(v=big$id[-1], lty=3, col="blue") abline(a=0, b=1, lty=3, col="blue") dev.off() #---------------------------------------------------------- # Residual Process Cusum Plot pdf(file="phuket3.pdf", height=8, width=8) par(mar=c(4, 4, 4.5, 0.2)) n <- nrow(x1$data) plot(residuals(x1) - 1:n, xlab="Event Number", ylab="Cusum") points(residuals(x0) - 1:n, lty=2, type="l", col="red") axis(3, at=big$id[-1], labels=big$date[-1], las=2) abline(v=big$id[-1], lty=3, col="blue") abline(a=0, b=0, lty=3, col="blue") dev.off() #---------------------------------------------------------- # Cusum Magnitude Plot pdf(file="phuket5.pdf", height=8, width=8) par(mar=c(4, 4, 4.5, 0.2)) lambda <- etas_gif(data=x1$data, evalpts=x1$data$time, params=x1$params[1:5]) plot(ts(cumsum(x1$data$magnitude-(1+sqrt(lambda)*x1$params[7])/x1$params[6])), ylab="Cusum", xlab="Event Number", ylim=c(-25, 20)) points(ts(cumsum(x0$data$magnitude-1/x0$params[6])), type="l", lty=2, col="red") abline(h=0, lty=3, col="blue") axis(3, at=big$id[-1], labels=big$date[-1], las=2) abline(v=big$id[-1], lty=3, col="blue") dev.off() #---------------------------------------------------------- # Simulation Example pdf(file="phuket4.pdf", height=8, width=8) par(mar=c(4, 4, 2, 0.2)) hist(y-x2$TT[1], breaks=seq(0, max(y-x2$TT[1]), length.out=20), xlab="Days Since 1 January 2009", main="") z <- quantile(y-x2$TT[1], probs=c(0.99, 0.95, 0.90, 0.80, 0.50)) abline(v=z, lty=2, col="blue") axis(3, at=z, labels=c(0.99, 0.95, 0.90, 0.80, 0.50)) box() dev.off() #---------------------------------------------------------- # Draw Map with Events data(Phuket) # magnitudes are rounded to 1 dp Phuket$magnitude <- Phuket$magnitude - 4.95 # Required graphics parameters usr <- c(89, 105, -5, 16) mai <- c(0.85, 0.85, 0.05, 0.15) # Determine required height of PDF file # given that we want a width of 6.5in width <- 6.5 aspect <- (usr[4]-usr[3])/(usr[2]-usr[1])/ cos(((usr[3]+usr[4])/2*pi)/180) height <- (width-mai[2]-mai[4])*aspect + mai[1] + mai[3] pdf(file="phuket6.pdf", height=height, width=width) # Layout graph boundaries and axes plot.new() par(mai=mai, usr=usr) title(xlab="Longitude") title(ylab="Latitude") axis(1) axis(2) box() # Add map to plot map("world2Hires", add=TRUE, interior=FALSE, col="gray35") # Add earthquake epicentres # Colour if within 30 days of one of 5 major events, else black for (j in 1:length(Phuket$latitude)) { if ((Phuket$time[j] > 360.04) & (Phuket$time[j] < 390.04)) col <- "red" else if ((Phuket$time[j] > 452.673) & (Phuket$time[j] < 482.673)) col <- "green" else if ((Phuket$time[j] > 1350.465) & (Phuket$time[j] < 1380.465)) col <- "blue" else if ((Phuket$time[j] > 1516.358) & (Phuket$time[j] < 1546.358)) col <- "magenta" else if ((Phuket$time[j] > 570.654) & (Phuket$time[j] < 600.654)) col <- "cyan" else col <- "black" points(Phuket$longitude[j], Phuket$latitude[j], cex=Phuket$magnitude[j]^1.5, col=col, lwd=Phuket$magnitude[j], pch = 1) } dev.off()