https://github.com/cran/fields
Tip revision: 6769ffc81115fbf0bf7d9c566cf7ac81be0049dc authored by Doug Nychka on 25 July 2005, 00:00:00 UTC
version 3.04
version 3.04
Tip revision: 6769ffc
sreg.r
"sreg" <-
function(x, y, lam = NA, df = NA, offset = 0, wt = rep(1, length(x)), cost = 1,
nstep.cv = 80, find.diagA = TRUE, trmin = 2.01, trmax =
length(unique(x)) * 0.95, lammin = NA, lammax = NA,
verbose = FALSE, do.cv = TRUE, method = "GCV", rmse = NA, lambda = NA)
{
call <- match.call()
out <- list()
out$call <- match.call()
class(out) <- c("sreg")
out$xraw <- x
out$yraw <- y
out$wt.raw <- wt
out$cost <- cost
out$offset <- offset
out$method <- method
if(length(x) != length(y))
stop(" X and Y do not match")
N <- length(y)
out$N <- N
out$nt <- 2
#
# figure out if the GCV function should be minimized
# and that value of lambda used for the estimate
#
if(!is.na(lambda[1])) {
lam <- lambda
}
#
#
if(is.na(lam[1]) & is.na(df[1])) {
eval.cv <- TRUE
find.min <- TRUE
}
else {
find.min <- FALSE
eval.cv <- FALSE
}
#
## find duplicate rows of the x vector
## first find integer tags to indicate replications
rep.info <- cat.matrix(matrix(x, ncol = 1))
out$rep.info <- rep.info
if(verbose) {
cat("rep.info", fill = TRUE)
print(rep.info)
}
if(max(rep.info) == N) {
shat.rep <- NA
shat.pure.error <- NA
out$pure.ss <- 0
YM <- y
weightsM <- wt
xM <- x[!duplicated(rep.info)]
}
else {
##
## do a simple 1-way ANOVA to get the replication error
##
rep.info.aov <- fast.1way(rep.info, y, wt)
shat.pure.error <- sqrt(rep.info.aov$MSE)
shat.rep <- shat.pure.error
YM <- rep.info.aov$means
weightsM <- rep.info.aov$w.means
xM <- as.matrix(x[!duplicated(rep.info)])
out$pure.ss <- rep.info.aov$SSE
if(verbose) {
cat(" rep info", fill = TRUE)
print(rep.info.aov)
}
}
#
# over write the real x's and y's with the collapsed group means.
# these orignal values are called xraw yraw.
#
out$np <- length(xM)
out$y <- c(YM)
out$x <- c(xM)
out$wt <- weightsM
out$shat.rep <- shat.rep
out$shat.pure.error <- shat.pure.error
xgrid <- sort(xM)
out$trace <- NA
#
# find lambda's if df's are given
#
if(!is.na(df[1])) {
lam <- rep(0, length(df))
for(k in 1:length(df)) {
lam[k] <- sreg.df.to.lambda(df[k], xM, weightsM)
}
}
if(do.cv) {
a <- gcv.sreg(out, lambda = lam, cost = cost, offset = offset,
nstep.cv = nstep.cv, verbose = verbose, trmin = trmin,
trmax = trmax, find.min = find.min, method = method,
rmse = rmse)
# if the spline should be evaluated at the GCV solutionwipe out lam grid
# and just use GCV lambda.
out$gcv.grid <- a$gcv.grid
out$lambda.est <- a$lambda.est
#
# save GCV estimate is that is what is needed
if(eval.cv) {
lam <- a$lambda.est[method, "lambda"]
out$trace <- a$lambda.est[method, "trA"]
out$shat.GCV <- a$lambda.est[method, "shat"]
}
#
# save trace and estimate of sigma an the traces
if(!find.min) {
out$trace <- c(a$gcv.grid[, "trA"])
out$shat.GCV <- c(a$gcv.grid[1, "shat"])
}
}
b <- list()
# lam can either be a grid or just the GCV value
NL <- length(lam)
NG <- length(xgrid)
h <- log(lam)
residuals <- matrix(NA, ncol = NL, nrow = length(out$yraw))
job <- as.integer(c(0, 3, 0))
predicted <- matrix(NA, ncol = NL, nrow = NG)
if(find.diagA) {
diagA <- matrix(NA, ncol = NL, nrow = out$np)
# add switch to find diag of A.
job <- as.integer(c(3, 3, 0))
}
for(k in 1:NL) {
#
# call cubic spline, note lambda is passed in log scale.
# spline solution evaluated at xgrid
#
b <- .Fortran("css",
h = as.double(h[k]),
npoint = as.integer(out$np),
x = as.double(out$x),
y = as.double(out$y),
wt = as.double(1/sqrt(out$wt)),
sy = as.double(rep(0, out$np)),
trace = as.double(0),
diag = as.double(c(cost, offset, rep(0, (out$np - 2)))),
cv = as.double(0),
ngrid = as.integer(NG),
xg = as.double(xgrid),
yg = as.double(rep(0, NG)),
job = as.integer(job),
ideriv = as.integer(0),
ierr = as.integer(0), PACKAGE="fields")
#
if(verbose) {
print(c(b$lambda, b$trace))
}
if(find.diagA) {
diagA[, k] <- b$diag
}
residuals[, k] <- out$yraw - splint(out$x, b$sy, out$xraw)
predicted[, k] <- b$yg
}
out$call <- call
out$lambda <- lam
out$do.cv <- do.cv
out$residuals <- residuals
out$fitted.values <- out$yraw - residuals
out$predicted <- list(x = xgrid, y = predicted)
if(length(lambda[1]) == 1) {
out$eff.df <- out$trace[1]
}
if(find.diagA) {
out$diagA <- diagA
}
out
}