# R/profile.R In SpatialExtremes: Modelling Spatial Extremes

```profile.maxstab <- function(fitted, param, range, n = 10,
plot = TRUE, conf = 0.90,
method = "RJ", square = "chol", ...){

if (!(method %in% c("RJ", "CB", "none")))
stop("'method' must be one of 'RJ', 'CB' or 'none'")

if (!(square %in% c("chol", "svd")))
stop("'square' must be one of 'chol' or 'svd'")

if (is.null(fitted\$std.err)){
warning("Profile confidence intervals cannot be computed whithout standard errors.")
method <- "none"
}

param.names <- names(fitted\$fitted.values)
n.param <- length(param.names)
idx.prof <- which(param.names == param)
fixed.param <- fitted\$fixed
fixed.names <- names(fixed.param)

ihessian <- fitted\$ihessian
var.cov <- fitted\$var.cov

if (method == "CB"){
##Chandler and Bate approach:
ivar.cov <- solve(var.cov)
hessian <- solve(ihessian)

if (square == "svd"){
svd.hessian <- svd(hessian)
svd.ivar.cov <- svd(ivar.cov)
M <- svd.hessian\$u %*% diag(sqrt(svd.hessian\$d)) %*%
t(svd.hessian\$u)
Madj <- svd.ivar.cov\$u %*% diag(sqrt(svd.ivar.cov\$d)) %*%
t(svd.ivar.cov\$u)
}

else{
M <- chol(hessian)
}

colnames(C) <- rownames(C) <- colnames(ihessian)
}

if (method == "RJ"){
##Rotnitzky and Jewell approach:
hessian <- solve(ihessian[param, param])
var.cov <- var.cov[param, param]

Q <- var.cov %*% hessian
eigen.val <- eigen(Q)\$values
}

start <- as.list(fitted\$fitted.values)
start <- start[-idx.prof]

fixed.values <- seq(range[1], range[2], length = n)
optfun <- fitted\$lik.fun

if (method == "CB"){
par <- NULL
for (fixed.val in range){
if (length(fixed.names) > 0)
body(optfun) <- parse(text = paste("nplk(", paste("p[",1:(n.param-1),"]", collapse = ","),
",", paste(fixed.names, "=", as.numeric(fixed.param),
collapse = ","), ",", param, "=", fixed.val, ")"))

else
body(optfun) <- parse(text = paste("nplk(", paste("p[",1:(n.param-1),"]", collapse = ","),
",", param, "=", fixed.val, ")"))

par <- rbind(par, optim(start, optfun)\$par)
}

par <- cbind(range, par)
colnames(par) <- c(param, param.names[-idx.prof])
par <- par[,names(fitted\$fitted)]
par.adj <- t(fitted(fitted) + C %*% (t(par) - fitted(fitted)))

}

llik <- rep(NA, n)
par <- matrix(NA, ncol = n.param - 1, nrow = n)

for (i in 1:n){
fixed.val <- fixed.values[i]

if (length(fixed.names) > 0)
body(optfun) <- parse(text = paste("nplk(", paste("p[",1:(n.param-1),"]", collapse = ","),
",", paste(fixed.names, "=", as.numeric(fixed.param),
collapse = ","), ",", param, "=", fixed.val, ")"))

else
body(optfun) <- parse(text = paste("nplk(", paste("p[",1:(n.param-1),"]", collapse = ","),
",", param, "=", fixed.val, ")"))

if (optfun(unlist(start)) >= 1e15)
reltol <- 1e-10

else
reltol <- 1e-6

opt <- optim(start, optfun, control = list(reltol = reltol, maxit = 10000))
llik[i] <- -opt\$value
par[i,] <- opt\$par
}

llik[llik <= -1e15] <- NA
ans <- cbind(fixed.values, llik, par)
colnames(ans) <- c(param, "llik", param.names[-idx.prof])
ans <- ans[,c("llik", names(fitted\$fitted))]

if (plot){
llikmax <- fitted\$logLik

if (method == "none")
plot(ans[,param], ans[,"llik"], xlab = param, ylab = "log-likelihood", ...)

if (method == "RJ"){
b.conf <- llikmax - .5 * eigen.val * qchisq(conf, 1)
plot(ans[,param], ans[,"llik"], xlab = param, ylab = "log-likelihood", ...)
abline(h = llikmax)
abline(h = b.conf)
}

if (method == "CB"){
ans[,-1] <- t(fitted(fitted) + solve(C) %*% (t(ans[,-1]) - fitted(fitted)))
b.conf <- llikmax - .5 * qchisq(conf, 1)
plot(ans[,param], ans[,"llik"], xlab = param, ylab = "log-likelihood", ...)
abline(h = llikmax)
abline(h = b.conf)
}
}

return(ans)

}

profile2d.maxstab <- function(fitted, params, ranges, n = 10,
plot = TRUE, ...){
param.names <- names(fitted\$fitted.values)
n.param <- length(param.names)
idx.prof <- which(param.names %in% params)
fixed.param <- fitted\$fixed
fixed.names <- names(fixed.param)

start <- as.list(fitted\$fitted.values)
start <- start[-idx.prof]

fixed.values1 <- seq(ranges[1,1], ranges[1,2], length = n)
fixed.values2 <- seq(ranges[2,1], ranges[2,2], length = n)

llik <- matrix(NA, n, n)

optfun <- fitted\$lik.fun

for (i in 1:n){
for (j in 1:n){
fixed.val <- c(fixed.values1[i], fixed.values2[j])
##We need to modify the body of optfun for each fixed value
if (length(fixed.names) > 0)
body(optfun) <- parse(text = paste("nplk(", paste("p[",1:(n.param-2),"]", collapse = ","),
",", paste(fixed.names, "=", as.numeric(fixed.param),
collapse = ","), ",",
paste(params, "=", fixed.val, collapse = ","), ")"))

else
body(optfun) <- parse(text = paste("nplk(", paste("p[",1:(n.param-2),"]", collapse = ","),
",", paste(params, "=", fixed.val, collapse = ","), ")"))

llik[i,j] <- -optim(start, optfun)\$value
}
}

ans <- list(coord = cbind(fixed.values1, fixed.values2), llik = llik)
colnames(ans\$coord) <- c(params)

if (plot)
filled.contour(fixed.values1, fixed.values2, llik,
xlab = params[1], ylab = params[2], ...)

return(ans)

}
```

## Try the SpatialExtremes package in your browser

Any scripts or data that you put into this service are public.

SpatialExtremes documentation built on May 2, 2019, 5:45 p.m.