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

```rbpspline <- function(y, x, knots, degree, penalty = "gcv", ...){

if ((degree %% 2) == 0)
stop("''degree'' must be an odd number")

if (!is.numeric(penalty) && (penalty != "cv") && (penalty != "gcv"))
stop("'penalty' must be either a numeric value, either 'cv' or 'gcv'")

if (is.numeric(penalty))
cv.hat <- NA

if (penalty == "cv"){
cv.fit <- cv(y, x, knots, degree, plot = FALSE, ...)
penalty <- cv.fit\$penalty
cv.hat <- cv.fit\$cv
}

if (penalty == "gcv"){
gcv.fit <- gcv(y, x, knots, degree, plot = FALSE, ...)
penalty <- gcv.fit\$penalty
cv.hat <- gcv.fit\$gcv
}

if (is.null(dim(x))){
idx <- order(x)
x <- x[idx]
y <- y[idx]
}

dsgn.mat <- rb(x, degree = degree, knots = knots,
penalty = penalty)
pen.mat <- dsgn.mat\$pen.mat
dsgn.mat <- dsgn.mat\$dsgn.mat

##To compute the beta vector, we use the Demmler-Reinch
##orthogonalization
chol.mat <- chol(crossprod(dsgn.mat))
chol.mat.inv <- solve(chol.mat)
svd.mat <- svd(t(chol.mat.inv) %*% pen.mat %*% chol.mat.inv)

A <- dsgn.mat %*% chol.mat.inv %*% svd.mat\$u
b <- t(A) %*% y

fitted.values <- A %*% (b / (1 + penalty^degree * svd.mat\$d))
smooth.mat <- A %*% diag(1 / (1 + penalty^degree * svd.mat\$d)) %*% t(A)
beta <- chol.mat.inv %*% svd.mat\$u %*% (b / (1 + penalty^degree * svd.mat\$d))

df <- sum(diag(smooth.mat))
res.df <- length(y) - 2 * sum(diag(smooth.mat)) +
sum(diag(smooth.mat %*% t(smooth.mat)))
rank <- sum(eigen(smooth.mat, only.values = TRUE)\$values != 0)

fitted <- list(x = x, obs = y, penalty = penalty, knots = knots,
degree = degree, y = fitted.values,
dsgn.mat = dsgn.mat, pen.mat = pen.mat, cv = cv.hat,
smooth.mat = smooth.mat, df = df, res.df = res.df,
rank = rank, beta = beta, call = match.call())

class(fitted) <- "pspline2"
return(fitted)
}

cv <- function(y, x, knots, degree, plot = TRUE, n.points = 150, ...){

dsgn.mat <- rb(x, degree = degree, knots = knots, penalty = NULL)
pen.mat <- dsgn.mat\$pen.mat
dsgn.mat <- dsgn.mat\$dsgn.mat

##To compute the beta vector, we use the Demmler-Reinch
##orthogonalization
chol.mat <- chol(crossprod(dsgn.mat))
chol.mat.inv <- solve(chol.mat)
svd.mat <- svd(t(chol.mat.inv) %*% pen.mat %*% chol.mat.inv)

A <- dsgn.mat %*% chol.mat.inv %*% svd.mat\$u
b <- t(A) %*% y

obj.fun <- function(penalty){

if (penalty < 0)
return(1e6)

fitted.values <- A %*% (b / (1 + penalty^degree * svd.mat\$d))
smooth.mat <- A %*% diag(1 / (1 + penalty^degree * svd.mat\$d)) %*% t(A)

cv <- sum(((y - fitted.values) / (1 - diag(smooth.mat)))^2)
return(cv)
}

opt <- nlm(obj.fun, sd(y), ...)

if (opt\$code >= 3)
warning("Optimization may have not succeeded")

if (plot){
cv.val <- NULL
lambda.vals <- seq(0, 2 * opt\$estimate, length = n.points)
for (lambda in lambda.vals)
cv.val <- c(cv.val, obj.fun(lambda))

cv.val[cv.val == 1e6] <- NA

plot(lambda.vals, cv.val, xlab = expression(lambda), ylab = "CV",
type = "l")
}

return(list(cv = opt\$minimum, penalty = opt\$estimate, nlm.code = opt\$code))

}

gcv <- function(y, x, knots, degree, plot = TRUE, n.points = 150, ...){

dsgn.mat <- rb(x, degree = degree, knots = knots, penalty = NULL)
pen.mat <- dsgn.mat\$pen.mat
dsgn.mat <- dsgn.mat\$dsgn.mat

##To compute the beta vector, we use the Demmler-Reinch
##orthogonalization
chol.mat <- chol(crossprod(dsgn.mat))
chol.mat.inv <- solve(chol.mat)
svd.mat <- svd(t(chol.mat.inv) %*% pen.mat %*% chol.mat.inv)

A <- dsgn.mat %*% chol.mat.inv %*% svd.mat\$u
b <- t(A) %*% y
n <- length(y)

obj.fun <- function(penalty){

if (penalty < 0)
return(1e6)

fitted.values <- A %*% (b / (1 + penalty^degree * svd.mat\$d))
smooth.mat <- A %*% diag(1 / (1 + penalty^degree * svd.mat\$d)) %*% t(A)

gcv <- n^2 * sum((y - fitted.values)^2) / (n - sum(diag(smooth.mat)))^2
return(gcv)
}

opt <- nlm(obj.fun, sd(y), ...)

if (opt\$code >= 3)
warning("Optimization may have not succeeded")

if (plot){
gcv.val <- NULL
lambda.vals <- seq(0, 2 * opt\$estimate, length = n.points)
for (lambda in lambda.vals)
gcv.val <- c(gcv.val, obj.fun(lambda))

gcv.val[gcv.val == 1e6] <- NA

plot(lambda.vals, gcv.val, xlab = expression(lambda), ylab = "GCV",
type = "l")
}

return(list(gcv = opt\$minimum, penalty = opt\$estimate, nlm.code = opt\$code))

}

rb <- function(..., knots, degree, penalty){

data <- cbind(...)

if ((degree %% 2) == 0)
stop("''degree'' must be an odd number")

if (missing(penalty))
penalty <- NULL

X0 <- NULL
for (i in 1:((degree - 1) / 2))
X0 <- cbind(X0, data^i)

X0 <- cbind(1, X0)

##Define the number of ``purely parametric'' parameters to be
##estimated i.e. weights related to radial basis function are not
##taken into account
n.ppar <- ncol(X0)

X1 <- NULL

if (is.null(dim(knots))){
n.knots <- length(knots)

for (k in 1:n.knots)
X1 <- cbind(X1, as.matrix(dist(c(knots[k], data)))[-1,1])

}

else{
n.knots <- nrow(knots)
for (k in 1:n.knots)
X1 <- cbind(X1, as.matrix(dist(rbind(knots[k,], data)))[-1,1])
}

X1 <- X1^degree
dsgn.mat <- cbind(X0, X1)

pen.mat <- matrix(0, nrow = dim(dsgn.mat)[2], ncol = dim(dsgn.mat)[2])
K <- as.matrix(dist(knots, upper = TRUE, diag = TRUE))^degree
pen.mat[-(1:n.ppar),-(1:n.ppar)] <- t(sqrt(K)) %*% sqrt(K)

dsgn.mat <- as.matrix(dsgn.mat)
pen.mat <- as.matrix(pen.mat)

ans <- list(dsgn.mat = dsgn.mat, pen.mat = pen.mat, degree = degree,
penalty = penalty, knots = knots, data = data, n.ppar = n.ppar,
call = match.call())

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.