# R/Pspline.q In pspline: Penalized Smoothing Splines

```smooth.Pspline <- function(x, y, w=rep(1, length(x)), norder=2,
df=norder+2, spar=0, method=1)
{

#  Computes order NORDER polynomial smoothing spline:  the spline
#    is piecewise of degree 2*NORDER - 1 and the norm of the
#    derivative of order is penalized

#  calls Fortran function spl

#  Arguments:
#  X       ...  argument values
#  Y       ...  N by NVAR matrix of function values to be smoothed
#  W       ...  weights (default all one's)
#  NORDER  ...  order of smoothing spline (default 2)
#  SPAR    ...  penalty parameter (default 0)
#  DF      ...  effective degrees of freedom (trace(hatmatrix))
#  METHOD  ...  smoothing method:  1  ...  fixed value of SPAR (default)
#                                  2  ...  fixed value of DF
#                                  3  ...  SPAR optimizes GCV criterion
#                                  4  ...  SPAR optimizes  CV criterion

#  Returns:  An object of class "smooth.Pspline" containing:
#  NORDER  ...  order of smoothing spline
#  X       ...  argument values
#  YSMTH   ...  N by NVAR matrix of values of the smoothed functions
#  LEV     ...  array of N leverage values
#  GCV     ...  generalized cross-validation coefficient
#  CV      ...  cross-validation coefficient
#  DF      ...  final effective degrees of freedom
#  SPAR    ...  final smoothing parameter value
#  MY.CALL ...  calling statement

my.call <- match.call()

n <- length(x)
if (is.matrix(y)) nvar <- ncol(y) else {
nvar <- 1
y <- as.matrix(y)
}
if (length(w) == 1) w <- rep(w,n)
if (nrow(y) != n | length(w) != n) stop("Argument arrays of wrong length")
if (method != 1 & method != 2 & method != 3 & method != 4) stop(
"Wrong value for METHOD")
if (norder <= 1 | norder >= 19) stop("Wrong value for NORDER")

yhat     <- matrix(0,n,nvar)
nworksiz <- (n-norder)*(4*norder + 3) + n
work     <- rep(0,nworksiz)
lev      <- rep(0,n)
gcv      <- 0
cv       <- 0
dfmax    <- n
ier      <- 0
irerun   <- 0

result <- .Fortran("pspline",
as.integer(n),      as.integer(nvar),  as.integer(norder),
as.double(x),       as.double(w),
as.double(y),       as.double(yhat),   as.double(lev),
as.double(gcv),     as.double(cv),     as.double(df),
as.double(spar),    as.double(dfmax),
as.double(work),    as.integer(method),
as.integer(irerun), as.integer(ier),
PACKAGE="pspline")

ier <- result[[17]]
if (ier == 1) stop ("N < 2*NORDER + 1")
if (ier == 2) stop ("NORDER < 2 or NORDER > 10")
if (ier == 3) stop ("NVAR < 1")
if (ier == 4) stop ("SPAR < 0")
if (ier == 5) stop ("X not strictly increasing")
if (ier == 6) stop ("W contains nonpositive values")
if (ier < 0 ) stop ("Singularity error in solving equations")

ysmth  <- matrix(result[[7]],n,nvar)
lev    <- result[[8]]
gcv    <- result[[9]]
cv     <- result[[10]]
df     <- result[[11]]
spar   <- result[[12]]
object <- list(norder = norder, x = x, ysmth = ysmth,  lev = lev,
gcv = gcv, cv = cv,  df = df,
spar = spar, call = my.call)
class(object) <- "smooth.Pspline"
object

}

predict.smooth.Pspline <- function(object, xarg, nderiv = 0, ...) {

if(missing(xarg)) return(object[c("x", "ysmth")])

x      <- object\$x
ysmth  <- object\$ysmth
norder <- 2*object\$norder
n    <- length(x)
nvar <- ncol(ysmth)
narg <- length(xarg)

if (nderiv < 0 | nderiv >= norder) stop("Violation of NDERIV >= NORDER.")

dy    <- matrix(0,narg,nvar)
work  <- rep(0,(2*norder+2)*n + norder)
ier   <- 0

result <- .Fortran("splifit",
as.integer(n),    as.integer(narg),
as.integer(nvar), as.integer(norder), as.integer(nderiv),
as.double(x),     as.double(ysmth),
as.double(xarg),  as.double(dy),
as.double(work),  as.integer(ier),
PACKAGE="pspline")

ier <- result[[11]]
if (ier == 1) stop (paste("N = ",n," not valid."))
if (ier == 2) stop ("A problem with knots detected.")
if (ier == 3) stop ("Singular coefficient matrix detected.")
if (ier == 4) stop (paste("NDERIV = ",nderiv," not valid."))
if (ier == 5) stop (paste("NORDER = ",norder," not valid."))
if (ier == 6) stop ("X values not strictly increasing.")

dy  <- matrix(result[[9]],narg,nvar)
return(dy)
}

plot.smooth.Pspline <- function(x, ...) {
if (is.vector(x\$ysmth) | dim(x\$ysmth)[[2]] == 1)
plot (x\$x, x\$ysmth, ...) else
matplot (x\$x, x\$ysmth, ...)
}

lines.smooth.Pspline <- function(x, ...) {
if (is.vector(x\$ysmth) | dim(x\$ysmth)[[2]] == 1)
lines (x\$x, x\$ysmth, ...) else
matlines (x\$x, x\$ysmth, ...)
}

print.smooth.Pspline <- function(x, ...)
{
if(!is.null(cl <- x\$call)) {
cat("Call:\n")
dput(cl)
}
cat("\nSmoothing Parameter (Spar):", format(x\$spar), "\n")
cat("Equivalent Degrees of Freedom (Df):", format(x\$df), "\n")
cat("GCV Criterion:", format(x\$gcv), "\n")
cat("CV  Criterion:", format(x\$cv),  "\n")
invisible(x)
}
```

## Try the pspline package in your browser

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

pspline documentation built on May 2, 2019, 4:01 a.m.