http://lib.stat.cmu.edu/general/gamfit
function (x, y = NULL, w = NULL, df, spar = NULL, cv = FALSE,
all.knots = FALSE, nknots = NULL, keep.data = TRUE, df.offset = 0,
penalty = 1, control.spar = list(), tol = 1e-06 * IQR(x))
{
contr.sp <- list(low = -1.5, high = 1.5, tol = 1e-04, eps = 2e-08,
maxit = 500, trace = getOption("verbose"))
contr.sp[names(control.spar)] <- control.spar
if (!all(sapply(contr.sp[1:4], is.numeric)) || contr.sp$tol <
0 || contr.sp$eps <= 0 || contr.sp$maxit <= 0)
stop("invalid 'control.spar'")
xy <- xy.coords(x, y)
y <- xy$y
x <- xy$x
if (!all(is.finite(c(x, y))))
stop("missing or infinite values in inputs are not allowed")
n <- length(x)
w <- if (is.null(w))
rep(1, n)
else {
if (n != length(w))
stop("lengths of 'x' and 'w' must match")
if (any(w < 0))
stop("all weights should be non-negative")
if (all(w == 0))
stop("some weights should be positive")
(w * sum(w > 0))/sum(w)
}
if (!is.finite(tol) || tol <= 0)
stop("'tol' must be strictly positive and finite")
xx <- round((x - mean(x))/tol)
nd <- !duplicated(xx)
ux <- sort(x[nd])
uxx <- sort(xx[nd])
nx <- length(ux)
if (nx <= 3L)
stop("need at least four unique 'x' values")
if (nx == n) {
ox <- TRUE
tmp <- cbind(w, w * y, w * y^2)[order(x), ]
}
else {
ox <- match(xx, uxx)
tapply1 <- function(X, INDEX, FUN = NULL, ..., simplify = TRUE) {
sapply(X = unname(split(X, INDEX)), FUN = FUN, ...,
simplify = simplify, USE.NAMES = FALSE)
}
tmp <- matrix(unlist(tapply1(seq_len(n), ox, function(i,
y, w) c(sum(w[i]), sum(w[i] * y[i]), sum(w[i] * y[i]^2)),
y = y, w = w), use.names = FALSE), ncol = 3, byrow = TRUE)
}
wbar <- tmp[, 1L]
ybar <- tmp[, 2L]/ifelse(wbar > 0, wbar, 1)
yssw <- sum(tmp[, 3L] - wbar * ybar^2)
if (is.na(cv) && !missing(df))
stop("'cv' must not be NA when 'df' is specified")
CV <- !is.na(cv) && cv
if (CV && nx < n)
warning("crossvalidation with non-unique 'x' values seems doubtful")
r.ux <- ux[nx] - ux[1L]
xbar <- (ux - ux[1L])/r.ux
if (all.knots) {
if (!is.null(nknots))
warning("'all.knots' is TRUE; 'nknots' specification is disregarded")
nknots <- nx
}
else {
if (is.null(nknots))
nknots <- n.knots(nx)
else if (!is.numeric(nknots))
stop("'nknots' must be numeric (in {1,..,n})")
else if (nknots < 1)
stop("'nknots' must be at least 1")
else if (nknots > nx)
stop("cannot use more inner knots than unique 'x' values")
}
knot <- c(rep(xbar[1], 3), if (all.knots) xbar else xbar[seq.int(1,
nx, length.out = nknots)], rep(xbar[nx], 3))
nk <- nknots + 2L
ispar <- if (is.null(spar) || missing(spar)) {
if (contr.sp$trace)
-1L
else 0L
}
else 1L
spar <- if (ispar == 1L)
as.double(spar)
else double(1)
icrit <- if (is.na(cv))
0L
else if (cv)
2L
else 1L
dofoff <- df.offset
if (!missing(df)) {
if (df > 1 && df <= nx) {
if (!missing(cv))
warning("specified both 'df' and 'cv'; will disregard the latter")
icrit <- 3L
dofoff <- df
}
else warning("you must supply 1 < df <= n, n = #{unique x} = ",
nx)
}
iparms <- as.integer(c(icrit, ispar, contr.sp$maxit))
names(iparms) <- c("icrit", "ispar", "iter")
keep.stuff <- FALSE
ans.names <- c("coef", "ty", "lev", "spar", "parms", "crit",
"iparms", "ier", if (keep.stuff) "scratch")
fit <- .Fortran(C_qsbart, as.double(penalty), as.double(dofoff),
x = as.double(xbar), y = as.double(ybar), w = as.double(wbar),
ssw = as.double(yssw), as.integer(nx), as.double(knot),
as.integer(nk), coef = double(nk), ty = double(nx), lev = double(if (is.na(cv)) 1L else nx),
crit = double(1), iparms = iparms, spar = spar, parms = unlist(contr.sp[1:4]),
isetup = as.integer(0), scratch = double(17L * nk + 1L),
ld4 = 4L, ldnk = 1L, ier = integer(1), DUP = FALSE)[ans.names]
wbar <- tmp[, 1]
if (is.na(cv))
lev <- df <- NA
else {
lev <- fit$lev
df <- sum(lev)
if (is.na(df))
stop("NA lev[]; probably smoothing parameter 'spar' way too large!")
}
if (fit$ier > 0L) {
sml <- fit$spar < 0.5
wtxt <- paste("smoothing parameter value too", if (sml)
"small"
else "large")
if (sml) {
stop(wtxt)
}
else {
fit$ty <- rep(mean(y), nx)
df <- 1
warning(wtxt, "\nsetting df = 1 __use with care!__")
}
}
cv.crit <- if (is.na(cv))
NA
else if (cv) {
ww <- wbar
ww[!(ww > 0)] <- 1
weighted.mean(((y - fit$ty[ox])/(1 - (lev[ox] * w)/ww[ox]))^2,
w)
}
else weighted.mean((y - fit$ty[ox])^2, w)/(1 - (df.offset +
penalty * df)/n)^2
pen.crit <- sum(wbar * (ybar - fit$ty)^2)
fit.object <- list(knot = knot, nk = nk, min = ux[1L], range = r.ux,
coef = fit$coef)
class(fit.object) <- "smooth.spline.fit"
object <- list(x = ux, y = fit$ty, w = wbar, yin = ybar,
data = if (keep.data) list(x = x, y = y, w = w), lev = lev,
cv.crit = cv.crit, pen.crit = pen.crit, crit = fit$crit,
df = df, spar = fit$spar, lambda = unname(fit$parms["low"]),
iparms = fit$iparms, fit = fit.object, call = match.call())
class(object) <- "smooth.spline"
object
}
