https://github.com/cran/fields
Tip revision: c71fb7f6ffa323303affebf0e35a0070faa9c24d authored by Doug Nychka on 10 May 2004, 00:00:00 UTC
version 1.7.2
version 1.7.2
Tip revision: c71fb7f
qsreg.family.R
"plot.qsreg" <-
function (out, pch = "*", main = NA)
{
old.par <- par("mfrow", "oma")
on.exit(par(old.par))
set.panel( 2, 2, relax=TRUE)
plot(out$x, out$y, xlab = "X", ylab = "y", pch = pch)
orderx <- order(out$x)
temp<- out$fitted.values[,c(out$ind.cv, out$ind.cv.ps)]
matlines(out$x[orderx], temp[orderx,], lty = 1, col = c(1,2))
##
# residual plot
#
matplot(out$x,
qsreg.psi( out$residuals[,c(out$ind.cv, out$ind.cv.ps)], out$alpha, out$sc)
, col=c(1,2), pch="o", ylab="Pseudo residuals", xlab="X")
yline(0)
if (nrow(out$cv.grid) > 1) {
ind <- out$cv.grid[, 3] < 1e+19
out$cv.grid <- out$cv.grid[ind, ]
matplot(out$cv.grid[, 2], cbind(out$cv.grid[, 3],
out$cv.grid[,6]), xlab = "Effective number of parameters",
ylab = "Log CV Rho function ",
log = "y", type="l", col=c(1,2))
xline( out$cv.grid[out$ind.cv, 2], col=1)
xline( out$cv.grid[out$ind.cv.ps, 2], col=2)
title(" CV curves",
cex = 0.5)
}
bplot(
qsreg.psi( out$residuals[,c(out$ind.cv, out$ind.cv.ps)], out$alpha, out$sc),
labels=c("CV", "CV pseudo")
)
yline( 0, col=2)
if (is.na(main))
mtext(deparse(out$call), cex = 1.3, outer = TRUE, line = -2)
else mtext(main, cex = 1.3, outer = TRUE, line = -2)
}
"predict.qsreg" <-
function (out, x, derivative = 0, model = out$ind.cv.ps)
{
if (missing(x))
x <- out$x
c(splint(out$predicted$x, out$predicted$y[, model], x, derivative = derivative))
}
"qsreg" <-
function (x, y, lam = NA, maxit = 50, maxit.cv = 10, tol = 1e-07,
offset = 0, sc = sqrt(var(y)) * 1e-05, alpha = 0.5, wt = rep(1,
length(x)), cost = 1, nstep.cv = 80, hmin = NA, hmax = NA,
trmin = 2 * 1.05, trmax = 0.95 * length(unique(x)))
{
out <- list()
class(out) <- c("qsreg")
N <- length(y)
out$N <- N
xgrid <- sort(unique(x))
if (length(x) != length(y))
stop(" X and Y do not match")
if (!is.na(lam[1]))
hgrid<- log( lam)
# find lambda grid
else
{
if (is.na(hmin)) {
hmin <- 0
for (k in 1:25) {
b <- qsreg.trace(lam = as.double(exp(hmin)),
x = x, y = y, wt = wt, cost = cost, maxit = maxit,
tol = tol, sc = sc, alpha = alpha)
if (b > trmax) {
break
}
hmin <- hmin - 1
}
}
if (is.na(hmax)) {
hmax <- 0
for (k in 1:25) {
b <- qsreg.trace(lam = as.double(exp(hmax)),
x = x, y = y, wt = wt, cost = cost, maxit = maxit,
tol = tol, sc = sc, alpha = alpha)
if (b < trmin) {
break
}
hmax <- hmax + 1
}
}
h<- seq( hmin,hmax,,nstep.cv)
lam<- exp( h)
}
# now loop through values for lam ( really log lam)
b <- list()
NL <- length(lam)
NG <- length(xgrid)
h <- log(lam)
residuals <- matrix(NA, ncol = NL, nrow = N)
diagA <- residuals
cv.ps <- rep(0, NL)
trace.ps <- rep(0, NL)
cv <- rep(0, NL)
predicted <- matrix(NA, ncol = NL, nrow = NG)
trace <- rep(0, NL)
converge <- rep(0, NL)
wt.old<- wt
for (k in 1:NL) {
b <- .Fortran("rcss", h = as.double(h[k]), npoint = as.integer(N),
x = as.double(x), y = as.double(y), wt = as.double(wt.old),
sy = as.double(rep(0, N)), trace = as.double(0),
diag = as.double(rep(0, N)), cv = as.double(0), ngrid = as.integer(NG),
xg = as.double(xgrid), yg = as.double(rep(0, NG)),
job = as.integer(c(3, 3, 0)), ideriv = as.integer(0),
din = as.double(c(cost, offset, maxit, tol, sc, alpha)),
dout = as.double(rep(0, 4)), ierr = as.integer(0), PACKAGE="fields")
residuals[, k] <- y - b$sy
diagA[, k] <- b$diag
cv[k] <- b$dout[4]
trace[k] <- b$trace
predicted[, k] <- b$yg
converge[k] <- b$dout[1]
wt.old<- b$wt
}
# second loop to find approx CV residuals based on pseudo values
y.pseudo<- rep(NA, N)
residuals.cv<- matrix( NA, ncol=NL, nrow= length( x))
for (k in 1:NL) {
y.pseudo<- (sc)*qsreg.psi(residuals[, k], alpha=alpha, C=sc) + y-residuals[,k]
#
# call the robust spline but set the cutoff for the huber weight so big
# it is essentially a LS spline this helps to match the lambda for robust spline
# with a LS one.
#
b <- .Fortran("rcss", h = as.double(h[k]), npoint = as.integer(N),
x = as.double(x), y = as.double(y.pseudo), wt = as.double(wt),
sy = as.double(rep(0, N)), trace = as.double(0),
diag = as.double(rep(0, N)), cv = as.double(0), ngrid = as.integer(NG),
xg = as.double(xgrid), yg = as.double(rep(0, NG)),
job = as.integer(c(3, 3, 0)), ideriv = as.integer(0),
din = as.double(c(cost, offset, maxit, tol, sqrt( var( y))*10, alpha)),
dout = as.double(rep(0, 4)), ierr = as.integer(0), PACKAGE="fields")
#
# CV residuals based on pseudo-data)
# Use the linear approximation Y_k - f.cv_k = (Y_k- f_k)/( 1- A_kk)
# f.cv_k = f_k/( 1- A_kk) - ( A_kk)Y_k/( 1- A_kk)
#
# Note: we find f.cv based on pseudo data but then consider its deviation
# from the actual data
#
f.cv<- (b$sy/(1-b$diag)) - b$diag*y.pseudo/(1-b$diag)
trace.ps[k]<- b$trace
residuals.cv[,k]<- (y-f.cv)
cv.ps[k]<- mean( qsreg.rho( y- f.cv, alpha=alpha, C=sc))
}
#
#
#
cv.grid <- cbind(lam, trace, cv, converge, trace.ps, cv.ps)
dimnames(cv.grid) <- list(NULL, c("lambda", "trace", "CV",
"iterations", "trace.PS", "CV.PS"))
#
ind.cv<- ( 1: NL) [ cv==min(cv)]
ind.cv.ps<- ( 1: NL) [ cv.ps==min(cv.ps)]
out$call <- match.call()
out$x <- x
out$y <- y
out$predicted<- list( x= xgrid, y= predicted)
out$trace<- trace
out$residuals.cv<- residuals.cv
out$residuals <- residuals
out$fitted.values<- y - residuals
out$cv.grid <- cv.grid
out$diagA <- diagA
out$sc <- sc
out$alpha <- alpha
out$ind.cv<- ind.cv
out$ind.cv.ps<- ind.cv.ps
out
}
"qsreg.fit" <-
function (x, y, lam, maxit = 50, maxit.cv = 10, tol = 1e-04,
offset = 0, sc = sqrt(var(y)) * 1e-07, alpha = 0.5, wt = rep(1,
length(x)), cost = 1)
{
N <- length(y)
if (length(x) != length(y))
stop(" X and Y do not match")
h <- log(lam)
temp <- .Fortran("rcss", h = as.double(log(lam)), npoint = as.integer(N),
x = as.double(x), y = as.double(y), wt = as.double(wt),
sy = as.double(rep(0, N)), trace = as.double(0), diag = as.double(rep(0,
N)), cv = as.double(0), ngrid = as.integer(0), xg = as.double(0),
yg = as.double(0), job = as.integer(c(3, 0, 0)), ideriv = as.integer(0),
din = as.double(c(cost, offset, maxit, tol, sc, alpha)),
dout = as.double(rep(0, 4)), ierr = as.integer(0), PACKAGE="fields")$dout
return(temp)
}
"qsreg.psi" <-
function( r,alpha=.5,C=1){
temp<- rep( NA, length( r))
r<- r/C
temp<- r
ind<- r>1
temp[ ind] <- 2*alpha
ind<- r<1& r>0
temp[ ind] <- (2*alpha*r[ind])
ind<- r< -1
temp[ ind] <- -2*(1-alpha)
ind<- r>-1& r<0
temp[ ind] <- 2*(1-alpha)*r[ind]
temp
}
"qsreg.rho" <-
function( r,alpha=.5,C=1){
temp<- rep( NA, length( r))
ind<- r>C
temp[ ind] <- 2*alpha*r[ind] - alpha*C
ind<- r<C& r>0
temp[ ind] <- (alpha*r[ind]**2)/C
ind<- r< -C
temp[ ind] <- -2*(1-alpha)*r[ind] - (1-alpha)*C
ind<- r>-C& r<0
temp[ ind] <- ((1-alpha)*r[ind]**2)/C
temp}