\name{isoreg}
\title{Isotonic / Monotone Regression}
\alias{isoreg}
%\alias{fitted.isoreg}
%\alias{residuals.isoreg}
\description{
  Compute the isotonic (monotonely increasing nonparametric) least
  squares regression which is piecewise constant.
}
\usage{
isoreg(x, y = NULL)
}
\arguments{
  \item{x, y}{in \code{isoreg}, coordinate vectors of the regression
    points.  Alternatively a single \dQuote{plotting} structure can be
    specified: see \code{\link{xy.coords}}.
  }
  \item{\dots}{potentially further arguments passed to methods.}
}
\details{
  The algorithm determines the convex minorant \eqn{m(x)} of the
  \emph{cumulative} data (i.e., \code{cumsum(y)}) which is piecewise
  linear and the result is \eqn{m'(x)}, a step function with level
  changes at locations where the convex \eqn{m(x)} touches the
  cumulative data polygon and changes slope.\cr
  \code{\link[stepfun]{as.stepfun}()} in \pkg{stepfun}
  returns a \code{\link[stepfun]{stepfun}} object
  which can be more parsimonious.
}
\value{
  \code{isoreg()} returns an object of class \code{isoreg} which is
  basically a list with components
  \item{x}{original (constructed) abscissa values \code{x}.}
  \item{y}{corresponding y values.}
  \item{yf}{fitted values corresponding to \emph{ordered} x values.}
  \item{yc}{cumulative y values corresponding to \emph{ordered} x values.}
  \item{iKnots}{integer vector giving indices where the fitted curve jumps,
    i.e., where the convex minorant has kinks.}
  \item{isOrd}{logical indicating if original x values were ordered
    increasingly already.}
  \item{ord}{\code{if(!isOrd)}: integer permutation \code{\link{order}(x)} of
    \emph{original} \code{x}.}
  \item{call}{the \code{\link{call}} to \code{isoreg()} used.}
}
\note{
  The code should be improved to accept \emph{weights} additionally and
  solve the corresponding weighted least squares problem.\cr
  \dQuote{Patches are welcome!}
}
\references{
  Barlow, R. E., Bartholomew, D. J., Bremner, J. M., and Brunk, H. D. (1972)
  \emph{Statistical inference under order restrictions}; Wiley, London.

  Robertson, T., Wright,F. T. and Dykstra, R. L. (1988)
  \emph{Order Restricted Statistical Inference}; Wiley, New York.
}
%%\author{Original C code by Brian Ripley; all else: Martin Maechler}
\seealso{the plotting method \code{\link{plot.isoreg}} with more examples;
  \code{\link[MASS]{isoMDS}()} from the \pkg{MASS} package internally
  uses isotonic regression.
}
\examples{
(ir <- isoreg(c(1,0,4,3,3,5,4,2,0)))
plot(ir, plot.type = "row")

(ir3 <- isoreg(y3 <- c(1,0,4,3,3,5,4,2, 3)))# last "3", not "0"
library(stepfun)
(fi3 <- as.stepfun(ir3))
(ir4 <- isoreg(1:10, y4 <- c(5, 9, 1:2, 5:8, 3, 8)))
cat("R^2 =", formatC(sum(residuals(ir4)^2) / (9*var(y4)), dig=2),"\n")
}
\keyword{regression}
\keyword{smooth}
\eof
\name{ksmooth}
\alias{ksmooth}
\title{Kernel Regression Smoother}
\description{
  The Nadaraya-Watson kernel regression estimate.
}
\usage{
ksmooth(x, y, kernel = c("box", "normal"), bandwidth = 0.5,
        range.x = range(x), n.points = max(100, length(x)), x.points)
}
\arguments{
 \item{x}{input x values}
 \item{y}{input y values}
 \item{kernel}{the kernel to be used.}
 \item{bandwidth}{the bandwidth. The kernels are scaled so that their
   quartiles (viewed as probability densities) are at
   \eqn{\pm}{+/-} \code{0.25*bandwidth}.}
 \item{range.x}{the range of points to be covered in the output.}
 \item{n.points}{the number of points at which to evaluate the fit.}
 \item{x.points}{points at which to evaluate the smoothed fit. If
   missing, \code{n.points} are chosen uniformly to cover \code{range.x}.}
}
\value{
  A list with components
  \item{x}{values at which the smoothed fit is evaluated. Guaranteed to
    be in increasing order.}
  \item{y}{fitted values corresponding to \code{x}.}
}
\note{This function is implemented purely for compatibility with S,
  although it is nowhere near as slow as the S function. Better kernel
  smoothers are available in other packages.}
\examples{
data(cars)
with(cars, {
    plot(speed, dist)
    lines(ksmooth(speed, dist, "normal", bandwidth=2), col=2)
    lines(ksmooth(speed, dist, "normal", bandwidth=5), col=3)
})
}
\keyword{smooth}
\eof
\name{loess}
\title{Local Polynomial Regression Fitting}
\usage{
loess(formula, data, weights, subset, na.action, model = FALSE,
      span = 0.75, enp.target, degree = 2,
      parametric = FALSE, drop.square = FALSE, normalize = TRUE,
      family = c("gaussian", "symmetric"),
      method = c("loess", "model.frame"),
      control = loess.control(\dots), \dots)
}
\alias{loess}
%\alias{anova.loess}
%\alias{print.loess}
%\alias{summary.loess}
%\alias{print.summary.loess}
\arguments{
  \item{formula}{a formula specifying the response and one to four
    numeric predictors (best specified via an interaction, but can also
    be specified additively).}
  \item{data}{an optional data frame within which to look first for the
    response, predictors and weights.}
  \item{weights}{optional weights for each case.}
  \item{subset}{an optional specification of a subset of the data to be
    used.}
  \item{na.action}{the action to be taken with missing values in the
    response or predictors.  The default is given by
    \code{getOption("na.action")}.}
  \item{model}{should the model frame be returned?}
  \item{span}{the parameter \eqn{\alpha} which controls the degree of
    smoothing.}
  \item{enp.target}{an alternative way to specify \code{span}, as the
    approximate equivalent number of parameters to be used.}
  \item{degree}{the degree of the polynomials to be used, up to 2.}
  \item{parametric}{should any terms be fitted globally rather than
    locally?  Terms can be specified by name, number or as a logical
    vector of the same length as the number of predictors.}
  \item{drop.square}{for fits with more than one predictor and
    \code{degree=2}, should the quadratic term (and cross-terms) be
    dropped for particular predictors?  Terms are specified in the same
    way as for \code{parametric}.}
  \item{normalize}{should the predictors be normalized to a common scale
    if there is more than one?  The normalization used is to set the
    10\% trimmed standard deviation to one.  Set to false for spatial
    coordinate predictors and others know to be a common scale.}
  \item{family}{if \code{"gaussian"} fitting is by least-squares, and if
    \code{"symmetric"} a re-descending M estimator is used with Tukey's
    biweight function.}
  \item{method}{fit the model or just extract the model frame.}
  \item{control}{control parameters: see \code{\link{loess.control}}.}
  \item{\dots}{control parameters can also be supplied directly.}
}
\description{
  Fit a polynomial surface determined by one or more numerical
  predictors, using local fitting.
}
\details{
  Fitting is done locally.  That is, for the fit at point \eqn{x}, the
  fit is made using points in a neighbourhood of \eqn{x}, weighted by
  their distance from \eqn{x} (with differences in \sQuote{parametric}
  variables being ignored when computing the distance). The size of the
  neighbourhood is controlled by \eqn{\alpha} (set by \code{span} or
  \code{enp.target}).  For \eqn{\alpha < 1}, the neighbourhood includes
  proportion \eqn{\alpha} of the points, and these have tricubic weighting
  (proportional to \eqn{(1 - \mathrm{(dist/maxdist)}^3)^3}{(1 -
    (dist/maxdist)^3)^3}. For \eqn{\alpha > 1}, all points are used, with
  the \sQuote{maximum distance} assumed to be
  \eqn{\alpha^{1/p}}{alpha^1/p} times the actual maximum distance for
  \eqn{p} explanatory variables.

  For the default family, fitting is by (weighted) least squares. For
  \code{family="symmetric"} a few iterations of an M-estimation
  procedure with Tukey's biweight are used. Be aware that as the initial
  value is the least-squares fit, this need not be a very resistant fit.

  It can be important to tune the control list to achieve acceptable
  speed. See \code{\link{loess.control}} for details.
}
\value{
  An object of class \code{"loess"}.
}
\references{
  W.S. Cleveland, E. Grosse and W.M. Shyu (1992) Local regression
  models. Chapter 8 of \emph{Statistical Models in S} eds J.M. Chambers
  and T.J. Hastie, Wadsworth & Brooks/Cole.
}
\author{
  B.D. Ripley, based on the \code{cloess} package of Cleveland,
  Grosse and Shyu avaliable at \url{http://www.netlib.org/a/}.
}
\note{
  As this is based on the \code{cloess} package available at
  \code{netlib}, it is similar to but not identical to the \code{loess}
  function of S. In particular, conditioning is not implemented.

  The memory usage of this implementation of \code{loess} is roughly
  quadratic in the number of points, with 1000 points taking about 10Mb.
}
\seealso{
  \code{\link{loess.control}},
  \code{\link{predict.loess}}.

  \code{\link[base]{lowess}}, the ancestor of \code{loess} (with
  different defaults!).
}
\examples{
data(cars)
cars.lo <- loess(dist ~ speed, cars)
predict(cars.lo, data.frame(speed = seq(5, 30, 1)), se = TRUE)
# to allow extrapolation
cars.lo2 <- loess(dist ~ speed, cars,
  control = loess.control(surface = "direct"))
predict(cars.lo2, data.frame(speed = seq(5, 30, 1)), se = TRUE)
}
\keyword{smooth}
\keyword{loess}
\eof
\name{loess.control}
\title{Set Parameters for Loess}
\usage{
loess.control(surface = c("interpolate", "direct"),
              statistics = c("approximate", "exact"),
              trace.hat = c("exact", "approximate"),
              cell = 0.2, iterations = 4, \dots)
}
\alias{loess.control}
\arguments{
 \item{surface}{should be fitted surface be computed exactly or via
   interpolation from a kd tree?}
 \item{statistics}{should the statistics be computed exactly or
   approximately? Exact computation can be very slow.}
 \item{trace.hat}{should the trace of the smoother matrix be computed
   exactly or approximately? It is recommended to use the approximation
   for more than about 1000 data points.}
 \item{cell}{if interpolation is used this controls the accuracy of the
   approximation via the maximum number of points in a  cell in the kd
   tree. Cells with more than \code{floor(n*span*cell)} points are subdivided.}
 \item{iterations}{the number of iterations used in robust fitting.}
  \item{\dots}{further arguments which are ignored.}
}
\description{
  Set control parameters for \code{loess} fits.
}
\value{
  A list with components
  \item{surface}{}
  \item{statistics}{}
  \item{trace.hat}{}
  \item{cell}{}
  \item{iterations}{}
  with meanings as explained under \sQuote{Arguments}.
}
\seealso{\code{\link{loess}}}
\keyword{smooth}
\keyword{loess}
\eof
\name{plot.isoreg}
\alias{plot.isoreg}
%\alias{print.isoreg}
\title{Plot Method for isoreg  Objects}
\description{
  The \code{\link{plot}} method for
  \R objects of class \code{\link{isoreg}}.
}
\usage{
\method{plot}{isoreg}(x, plot.type = c("single", "row.wise", "col.wise"),
      main = paste("Isotonic regression", deparse(x$call)),
      main2 = "Cumulative Data and Convex Minorant",
      xlab = "x0", ylab = "x$y",
      par.fit = list(col = "red", cex = 1.5, pch = 13, lwd = 1.5),
      mar = if (both) 0.1 + c(3.5, 2.5, 1, 1) else par("mar"),
      mgp = if (both) c(1.6, 0.7, 0) else par("mgp"),
      grid = length(x$x) < 12, \dots)
%\method{print}{isoreg}(x, digits = getOption("digits"), \dots)
}
\arguments{
  \item{x}{an \code{\link{isoreg}} object.}
  \item{plot.type}{character indicating which type of plot is desired.
    The first (default) only draws the data and the fit, where the
    others add a plot of the cumulative data and fit.}
  \item{main}{main title of plot, see \code{\link{title}}.}
  \item{main2}{title for second (cumulative) plot.}
  \item{xlab, ylab}{x- and y- axis annotation.}
  \item{par.fit}{a \code{\link{list}} of arguments (for
    \code{\link{points}} and \code{\link{lines}}) for drawing the fit.}
  \item{mar, mgp}{graphical parameters, see \code{\link{par}}, mainly
    for the case of two plots.}
  \item{grid}{logical indicating if grid lines should be drawn.  If
    true, \code{\link{grid}()} is used for the first plot, where as
    vertical lines are drawn at \dQuote{touching} points for the
    cumulative plot.}
  \item{\dots}{further arguments passed to and from methods.}
%  \item{digits}{integer indicating number of digits to use for printing.}
}
\seealso{\code{\link{isoreg}} for computation of \code{isoreg} objects.}
\examples{
example(isoreg) # for the examples there

## 'same' plot as above, "proving" that only ranks of 'x' are important
plot(isoreg(2^(1:9), c(1,0,4,3,3,5,4,2,0)), plot.t = "row", log = "x")

plot(ir3, plot.type = "row", ylab = "y3")
plot(isoreg(y3 - 4), plot.t="r", ylab = "y3 - 4")
plot(ir4, plot.type = "ro",  ylab = "y4", xlab = "x = 1:n")

## experiment a bit with these (C-c C-j):
plot(isoreg(sample(9),  y3), plot.type="row")
plot(isoreg(sample(9),  y3), plot.type="col.wise")

plot(ir <- isoreg(sample(10), sample(10, replace = TRUE)), plot.t = "r")
}
\keyword{hplot}
\keyword{print}
\eof
% file modreg/man/plot.ppr.Rd
% copyright (C) 1995-8 B. D. Ripley
%
\name{plot.ppr}
\alias{plot.ppr}
\title{
Plot Ridge Functions for Projection Pursuit Regression Fit
}
\description{
Plot ridge functions for projection pursuit regression fit.
}
\usage{
\method{plot}{ppr}(x, ask, type = "o", \dots)
}
\arguments{
  \item{x}{
    A fit of class \code{"ppr"} as produced by a call to \code{ppr}.
  }
  \item{ask}{
    the graphics parameter \code{ask}: see \code{par} for details.
    If set to \code{TRUE} will ask between the plot of each cross-section.
  }
  \item{type}{
    the type of line to draw
  }
  \item{\dots}{
    further graphical parameters
  }
}
\value{
  None
}
\section{Side Effects}{
  A series of plots are drawn on the current graphical device, one for
  each term in the fit.
}
\seealso{
  \code{\link{ppr}}, \code{\link{par}}
}
\examples{
data(rock)
attach(rock)
area1 <- area/10000; peri1 <- peri/10000
par(mfrow=c(3,2))# maybe: , pty="s")
rock.ppr <- ppr(log(perm) ~ area1 + peri1 + shape,
                data = rock, nterms = 2, max.terms = 5)
plot(rock.ppr, main="ppr(log(perm)~ ., nterms=2, max.terms=5)")
plot(update(rock.ppr, bass=5), main = "update(..., bass = 5)")
plot(update(rock.ppr, sm.method="gcv", gcvpen=2),
     main = "update(..., sm.method=\"gcv\", gcvpen=2)")
detach()
}
\keyword{hplot}
% Converted by Sd2Rd version 0.2-a5.
\eof
% file modreg/man/ppr.Rd
% copyright (C) 1995-8 B. D. Ripley
% copyright (C) 2000   The R Core Development Team
%
\name{ppr}
\alias{ppr}
\alias{ppr.default}
\alias{ppr.formula}
%\alias{predict.ppr}
%\alias{print.ppr}
%\alias{summary.ppr}
%\alias{print.summary.ppr}
\title{Projection Pursuit Regression}
\description{
  Fit a projection pursuit regression model.
}
\usage{
ppr(x, \dots)

\method{ppr}{formula}(formula, data = sys.parent(), weights,
    subset, na.action, contrasts = NULL, \dots)

\method{ppr}{default}(x, y, weights = rep(1,n),
    ww = rep(1,q), nterms, max.terms = nterms, optlevel = 2,
    sm.method = c("supsmu", "spline", "gcvspline"),
    bass = 0, span = 0, df = 5, gcvpen = 1, \dots)
}
\arguments{
  \item{formula}{
    a formula specifying one or more response variables and the
    explanatory variables.
  }
  \item{x}{
    matrix of explanatory variables.  Rows represent observations, and
    columns represent variables.  Missing values are not accepted.
  }
  \item{y}{
    matrix of response variables.  Rows represent observations, and
    columns represent variables.  Missing values are not accepted.
  }
  \item{nterms}{number of terms to include in the final model.}
  \item{data}{
    data frame from which variables specified in \code{formula} are
    preferentially to be taken.
  }
  \item{weights}{a vector of weights \code{w_i} for each \emph{case}.}
  \item{ww}{
    a vector of weights for each \emph{response}, so the fit criterion is
    the sum over case \code{i} and responses \code{j} of
    \code{w_i ww_j (y_ij - fit_ij)^2} divided by the sum of \code{w_i}.
  }
  \item{subset}{
    an index vector specifying the cases to be used in the training
    sample.  (NOTE: If given, this argument must be named.)
  }
  \item{na.action}{
    a function to specify the action to be taken if \code{\link{NA}}s are
    found. The default action is given by \code{getOption("na.action")}.
    (NOTE: If given, this argument must be named.)
  }
  \item{contrasts}{
    the contrasts to be used when any factor explanatory variables are coded.
  }
  \item{max.terms}{
    maximum number of terms to choose from when building the model.
  }
  \item{optlevel}{
    integer from 0 to 3 which determines the thoroughness of an
    optimization routine in the SMART program. See the \bold{Details} section.
  }
  \item{sm.method}{
    the method used for smoothing the ridge functions.  The default is to
    use Friedman's super smoother \code{supsmu}.  The alternatives are to use
    the smoothing spline code underlying \code{smooth.spline}, either with a
    specified (equivalent) degrees of freedom for each ridge functions, or
    to allow the smoothness to be chosen by GCV.
  }
  \item{bass}{
    super smoother bass tone control used with automatic span selection
    (see \code{supsmu}); the range of values is 0 to 10, with larger values
    resulting in increased smoothing.
  }
  \item{span}{
    super smoother span control (see \code{supsmu}).  The default, \code{0},
    results in automatic span selection by local cross validation. \code{span}
    can also take a value in \code{(0, 1]}.
  }
  \item{df}{
    if \code{sm.method} is \code{"spline"} specifies the smoothness of
    each ridge term via the requested equivalent degrees of freedom.
  }
  \item{gcvpen}{
    if \code{sm.method} is \code{"gcvspline"} this is the penalty used
    in the GCV selection for each degree of freedom used.
  }
  \item{\dots}{arguments to be passed to or from other methods.}
}
\value{
A list with the following components, many of which are for use by the
method functions.

\item{call}{the matched call}
\item{p}{the number of explanatory variables (after any coding)}
\item{q}{the number of response variables}
\item{mu}{the argument \code{nterms}}
\item{ml}{the argument \code{max.terms}}
\item{gof}{the overall residual (weighted) sum of squares for the
  selected model}
\item{gofn}{the overall residual (weighted) sum of squares against the
  number of terms, up to \code{max.terms}.  Will be invalid (and zero)
  for less than \code{nterms}.}
\item{df}{the argument \code{df}}
\item{edf}{if \code{sm.method} is \code{"spline"} or \code{"gcvspline"}
  the equivalent number of degrees of freedom for each ridge term used.}
\item{xnames}{the names of the explanatory variables}
\item{ynames}{the names of the response variables}
\item{alpha}{a matrix of the projection directions, with a column for
  each ridge term}
\item{beta}{a matrix of the coefficients applied for each response to
  the ridge terms: the rows are the responses and the columns the ridge terms}
\item{yb}{the weighted means of each response}
\item{ys}{the overall scale factor used: internally the responses are
  divided by \code{ys} to have unit total weighted sum of squares.}
\item{fitted.values}{the fitted values, as a matrix if \code{q > 1}.}
\item{residuals}{the residuals, as a matrix if \code{q > 1}.}
\item{smod}{internal work array, which includes the ridge functions
  evaluated at the training set points.}
}
\details{
  The basic method is given by Friedman (1984), and is essentially the
  same code used by S-PLUS's \code{ppreg}.  This code is extremely
  sensitive to the compiler used.

  The algorithm first adds up to \code{max.terms} ridge terms one at a
  time; it will use less if it is unable to find a term to add that makes
  sufficient difference.  It then removes the least
  \dQuote{\emph{important}} term at each step until \code{nterm} terms
  are left.

  The levels of optimization (argument \code{optlevel})
  differ in how thoroughly the models are refitted during this process.
  At level 0 the existing ridge terms are not refitted.  At level 1
  the projection directions are not refitted, but the ridge
  functions and the regression coefficients are.
%
  Levels 2 and 3 refit all the terms and are equivalent for one
  response; level 3 is more careful to re-balance the contributions
  from each regressor at each step and so is a little less likely to
  converge to a saddle point of the sum of squares criterion.
}
\references{
  Friedman, J. H. and Stuetzle, W. (1981)
  Projection pursuit regression.
  \emph{Journal of the American Statistical Association},
  \bold{76}, 817--823.

  Friedman, J. H. (1984)
  SMART User's Guide.
  Laboratory for Computational Statistics, Stanford University Technical
  Report No. 1.

  Venables, W. N. & Ripley, B. D. (2002)
  \emph{Modern Applied Statistics with S.}  Springer.
}
\seealso{
  \code{\link{plot.ppr}}, \code{\link{supsmu}}, \code{\link{smooth.spline}}
}
\examples{
# Note: your numerical values may differ
data(rock)
attach(rock)
area1 <- area/10000; peri1 <- peri/10000
rock.ppr <- ppr(log(perm) ~ area1 + peri1 + shape,
                data = rock, nterms = 2, max.terms = 5)
rock.ppr
# Call:
# ppr.formula(formula = log(perm) ~ area1 + peri1 + shape, data = rock,
#     nterms = 2, max.terms = 5)
#
# Goodness of fit:
#  2 terms  3 terms  4 terms  5 terms
# 8.737806 5.289517 4.745799 4.490378

summary(rock.ppr)
# .....  (same as above)
# .....
#
# Projection direction vectors:
#       term 1      term 2
# area1  0.34357179  0.37071027
# peri1 -0.93781471 -0.61923542
# shape  0.04961846  0.69218595
#
# Coefficients of ridge terms:
#    term 1    term 2
# 1.6079271 0.5460971

par(mfrow=c(3,2))# maybe: , pty="s")
plot(rock.ppr, main="ppr(log(perm)~ ., nterms=2, max.terms=5)")
plot(update(rock.ppr, bass=5), main = "update(..., bass = 5)")
plot(update(rock.ppr, sm.method="gcv", gcvpen=2),
     main = "update(..., sm.method=\"gcv\", gcvpen=2)")
cbind(perm=rock$perm, prediction=round(exp(predict(rock.ppr)), 1))
detach()
}
\keyword{regression}

\eof
\name{predict.loess}
\alias{predict.loess}
\title{Predict Loess Curve or Surface}
\description{
 Predictions from a \code{loess} fit, optionally with standard errors.
}
\usage{
\method{predict}{loess}(object, newdata = NULL, se = FALSE, \dots)
}
\arguments{
 \item{object}{an object fitted by \code{loess}.}
 \item{newdata}{an optional data frame specifying points at which to do
   the predictions. If missing, the original data points are used.}
 \item{se}{should standard errors be computed?}
 \item{\dots}{arguments passed to or from other methods.}
}
\details{
  The standard errors calculation is slower than prediction.

  When the fit was made using \code{surface="interpolate"} (the
  default), \code{predict.loess} will not extrapolate -- so points outside
  an axis-aligned hypercube enclosing the original data will have
  missing (\code{NA}) predictions and standard errors.
}
\value{
  If \code{se = FALSE}, a vector giving the prediction for each row of
  \code{newdata} (or the original data). If \code{se = TRUE}, a list
  containing components
  \item{fit}{the predicted values.}
  \item{se}{an estimated standard error for each predicted value.}
  \item{residual.scale}{the estimated scale of the residuals used in
    computing the standard errors.}
  \item{df}{an estimate of the effective degrees of freedom used in
    estimating the residual scale, intended for use with t-based
    confidence intervals. }
}
\author{
  B.D. Ripley, based on the \code{cloess} package of Cleveland,
  Grosse and Shyu.
}

\seealso{\code{\link{loess}}}

\examples{
data(cars)
cars.lo <- loess(dist ~ speed, cars)
predict(cars.lo, data.frame(speed=seq(5, 30, 1)), se=TRUE)
# to get extrapolation
cars.lo2 <- loess(dist ~ speed, cars,
  control=loess.control(surface="direct"))
predict(cars.lo2, data.frame(speed=seq(5, 30, 1)), se=TRUE)
}
\keyword{smooth}
\eof
\name{predict.smooth.spline}
\alias{predict.smooth.spline}
%\alias{predict.smooth.spline.fit}
\title{Predict from Smoothing Spline Fit}
\description{
  Predict a smoothing spline fit at new points, return the derivative if
  desired. The predicted fit is linear beyond the original data.
}
\usage{
\method{predict}{smooth.spline}(object, x, deriv = 0, \dots)
}
\arguments{
  \item{object}{a fit from \code{smooth.spline}.}
  \item{x}{the new values of x.}
  \item{deriv}{integer; the order of the derivative required.}
  \item{\dots}{further arguments passed to or from other methods.}
}
\value{
  A list with components
  \item{x}{The input \code{x}.}
  \item{y}{The fitted values or derivatives at \code{x}.}
}
\seealso{\code{\link{smooth.spline}}
}
\examples{
data(cars)
attach(cars)
cars.spl <- smooth.spline(speed, dist, df=6.4)
\dontshow{
print.default(cars.spl)
}

## "Proof" that the derivatives are okay, by comparing with approximation
diff.quot <- function(x,y) {
  ## Difference quotient (central differences where available)
  n <- length(x); i1 <- 1:2; i2 <- (n-1):n
  c(diff(y[i1]) / diff(x[i1]), (y[-i1] - y[-i2]) / (x[-i1] - x[-i2]),
    diff(y[i2]) / diff(x[i2]))
}

xx <- unique(sort(c(seq(0,30, by = .2), kn <- unique(speed))))
i.kn <- match(kn, xx)# indices of knots within xx
op <- par(mfrow = c(2,2))
plot(speed, dist, xlim = range(xx), main = "Smooth.spline & derivatives")
lines(pp <- predict(cars.spl, xx), col = "red")
points(kn, pp$y[i.kn], pch = 3, col="dark red")
mtext("s(x)", col = "red")
for(d in 1:3){
  n <- length(pp$x)
  plot(pp$x, diff.quot(pp$x,pp$y), type = 'l', xlab="x", ylab="",
       col = "blue", col.main = "red",
       main= paste("s",paste(rep("'",d), collapse=""),"(x)", sep=""))
  mtext("Difference quotient approx.(last)", col = "blue")
  lines(pp <- predict(cars.spl, xx, deriv = d), col = "red")
\dontshow{
  print(pp)
}
  points(kn, pp$y[i.kn], pch = 3, col="dark red")
  abline(h=0, lty = 3, col = "gray")
}
detach(); par(op)
}
\keyword{smooth}
\eof
% file MASS/rock.d
% copyright (C) 1994-8 W. N. Venables and B. D. Ripley
%
% Manually converted to Rd data file description format by KH.
\name{rock}
\docType{data}
\alias{rock}
\title{Measurements on Petroleum Rock Samples}
\description{
    Measurements on 48 rock samples from a petroleum reservoir.}
\usage{data(rock)}
\format{
    A data frame with 48 rows and 4 numeric columns.

    \tabular{rll}{
	[,1] \tab area  \tab area of pores space, in pixels
	out of 256 by 256 \cr
	[,2] \tab peri  \tab perimeter in pixels \cr
	[,3] \tab shape \tab perimeter/sqrt(area) \cr
	[,4] \tab perm  \tab permeability in milli-Darcies
    }
}
\details{
    Twelve core samples from petroleum reservoirs were sampled by 4
    cross-sections.  Each core sample was measured for permeability, and
    each cross-section has total area of pores, total perimeter of
    pores, and shape.
}
\source{
    Data from BP Research, image analysis by Ronit Katz, U. Oxford.
}
\keyword{datasets}
\eof
\name{runmed}
\title{Running Medians -- Robust Scatter Plot Smoothing}
\alias{runmed}
\description{
  Compute running medians of odd span.  This is the \dQuote{most robust}
  scatter plot smoothing possible.  For efficiency (and historical
  reason), you can use one of two different algorithms giving identical
  results.
}
\usage{
runmed(x, k, endrule = c("median","keep","constant"),
       algorithm = NULL, print.level = 0)
}
\arguments{
  \item{x}{numeric vector, the \dQuote{dependent} variable to be
    smoothed.}
  \item{k}{integer width of median window; must be odd.  Turlach had a
    default of \code{k <- 1 + 2 * min((n-1)\%/\% 2, ceiling(0.1*n))}.
    Use \code{k = 3} for \dQuote{minimal} robust smoothing eliminating
    isolated outliers.}
  \item{endrule}{character string indicating how the values at the
    beginning and the end (of the data) should be treated.
    \describe{
      \item{\code{"keep"}}{keeps the first and last \eqn{k_2}{k2} values
	at both ends, where \eqn{k_2}{k2} is the half-bandwidth \code{k2
	  = k \%/\% 2},
	i.e., \code{y[j] = x[j]} for \eqn{j \in \{1,\ldots,k_2;
	  n-k_2+1,\ldots,n\}}{j = 1,..,k2 and (n-k2+1),..,n};}
      \item{\code{"constant"}}{copies \code{median(y[1:k2])} to the first
	values and analogously for the last ones making the smoothed ends
	\emph{constant};}
      \item{\code{"median"}}{the default, smoothes the ends by using
	symmetrical medians of subsequently smaller bandwidth, but for
	the very first and last value where Tukey's robust end-point
	rule is applied, see \code{\link{smoothEnds}}.}
    }
  }
  \item{algorithm}{character string (partially matching \code{"Turlach"} or
    \code{"Stuetzle"}) or the default \code{NULL}, specifying which algorithm
    should be applied.  The default choice depends on \code{n = length(x)}
    and \code{k} where \code{"Turlach"} will be used for larger problems.}
  \item{print.level}{integer, indicating verboseness of algorithm;
    should rarely be changed by average users.}
}
\value{vector of smoothed values of the same length as \code{x} with an
  \code{\link{attr}}ibute \code{k} containing (the \sQuote{oddified})
  \code{k}.
}
\details{
  Apart from the end values, the result \code{y = runmed(x, k)} simply has
  \code{y[j] = median(x[(j-k2):(j+k2)])} (k = 2*k2+1), computed very
  efficiently.

  The two algorithms are internally entirely different:
  \describe{
    \item{"Turlach"}{is the Härdle-Steiger
      algorithm (see Ref.) as implemented by Berwin Turlach.
      A tree algorithm is used, ensuring performance \eqn{O(n \log
	k)}{O(n * log(k))} where \code{n <- length(x)} which is
      asymptotically optimal.}
    \item{"Stuetzle"}{is the (older) Stuetzle-Friedman implementation
      which makes use of median \emph{updating} when one observation
      enters and one leaves the smoothing window.  While this performs as
      \eqn{O(n \times k)}{O(n * k)} which is slower asymptotically, it is
      considerably faster for small \eqn{k} or \eqn{n}.}
  }
}
\references{
  Härdle, W. and Steiger, W. (1995)
  [Algorithm AS 296] Optimal median smoothing,
  \emph{Applied Statistics} \bold{44}, 258--264.

  Jerome H. Friedman and Werner Stuetzle (1982)
  \emph{Smoothing of Scatterplots};
  Report, Dep. Statistics, Stanford U., Project Orion 003.

  Martin Maechler (2003)
  Fast Running Medians: Finite Sample and Asymptotic Optimality;
  working paper available from the author.
}
\author{Martin Maechler \email{maechler@stat.math.ethz.ch},
  based on Fortran code from Werner Stuetzle and S-plus and C code from
  Berwin Turlach.
}
\seealso{
  \code{\link{smoothEnds}} which implements Tukey's end point rule and
  is called by default from \code{runmed(*, endrule = "median")}.
  \code{\link[eda]{smooth}} (from \pkg{eda} package) uses running
  medians of 3 for its compound smoothers.
}
\examples{
example(nhtemp)#> data(nhtemp)
myNHT <- as.vector(nhtemp)
myNHT[20] <- 2 * nhtemp[20]
plot(myNHT, type="b", ylim = c(48,60), main = "Running Medians Example")
lines(runmed(myNHT, 7), col = "red")

## special: multiple y values for one x
data(cars)
plot(cars, main = "'cars' data and runmed(dist, 3)")
lines(cars, col = "light gray", type = "c")
with(cars, lines(speed, runmed(dist, k = 3), col = 2))
%% FIXME: Show how to do it properly ! tapply(*, unique(.), median)

## nice quadratic with a few outliers
y <- ys <- (-20:20)^2
y [c(1,10,21,41)] <- c(150, 30, 400, 450)
all(y == runmed(y, 1)) # 1-neigborhood <==> interpolation
plot(y) ## lines(y, lwd=.1, col="light gray")
lines(lowess(seq(y),y, f = .3), col = "brown")
lines(runmed(y, 7), lwd=2, col = "blue")
lines(runmed(y,11), lwd=2, col = "red")

## Lowess is not robust
y <- ys ; y[21] <- 6666 ; x <- seq(y)
col <- c("black", "brown","blue")
plot(y, col=col[1])
lines(lowess(x,y, f = .3), col = col[2])
%% predict(loess(y ~ x, span = 0.3, degree=1, family = "symmetric"))
%% gives 6-line warning but does NOT break down
lines(runmed(y, 7),      lwd=2, col = col[3])
legend(length(y),max(y), c("data", "lowess(y, f = 0.3)", "runmed(y, 7)"),
       xjust = 1, col = col, lty = c(0, 1,1), pch = c(1,NA,NA))
}
\keyword{smooth}
\keyword{robust}
\eof
\name{scatter.smooth}
\title{Scatter Plot with Smooth Curve Fitted by Loess}
\alias{scatter.smooth}
\alias{loess.smooth}
\description{
  Plot and add a smooth curve computed by \code{loess} to a scatter plot.
}
\usage{
scatter.smooth(x, y, span = 2/3, degree = 1,
    family = c("symmetric", "gaussian"),
    xlab = deparse(substitute(x)), ylab = deparse(substitute(y)),
    ylim = range(y, prediction$y), evaluation = 50, \dots)

loess.smooth(x, y, span = 2/3, degree = 1,
    family = c("symmetric", "gaussian"), evaluation = 50, \dots)
}
\arguments{
  \item{x}{x coordinates for scatter plot.}
  \item{y}{y coordinates for scatter plot.}
  \item{span}{smoothness parameter for \code{loess}.}
  \item{degree}{degree of local polynomial used.}
  \item{family}{if \code{"gaussian"} fitting is by least-squares, and if
    \code{family="symmetric"} a re-descending M estimator is used.} 
  \item{xlab}{label for x axis.}
  \item{ylab}{label for y axis.}
  \item{ylim}{the y limits of the plot.}
  \item{evaluation}{number of points at which to evaluate the smooth
    curve.}
  \item{\dots}{graphical parameters.}
}
\details{
  \code{loess.smooth} is an auxiliary function which evaluates the
  \code{loess} smooth at \code{evaluation} equally spaced points
  covering the range of \code{x}.
}
\value{
  For \code{scatter.smooth}, none.

  For \code{loess.smooth}, a list with two components, \code{x} (the
  grid of evaluation points) and \code{y} (the smoothed values at the
  grid points).
}
\seealso{\code{\link{loess}}}
\examples{
data(cars)
attach(cars)
scatter.smooth(speed, dist)
detach()
}
\keyword{smooth}
\eof
\name{smooth.spline}
\alias{smooth.spline}
%\alias{print.smooth.spline}
\title{Fit a Smoothing Spline}
\description{
  Fits a cubic smoothing spline to the supplied data.
}
\usage{
smooth.spline(x, y = NULL, w = NULL, df, spar = NULL,
              cv = FALSE, all.knots = FALSE, nknots = NULL,
              df.offset = 0, penalty = 1, control.spar = list())
}
\arguments{
 \item{x}{a vector giving the values of the predictor variable, or  a
  list or a two-column matrix specifying x and y. }
 \item{y}{responses. If \code{y} is missing, the responses are assumed
   to be specified by \code{x}.}
 \item{w}{optional vector of weights of the same length as \code{x};
   defaults to all 1.}
 \item{df}{the desired equivalent number of degrees of freedom (trace of
   the smoother matrix).}
 \item{spar}{smoothing parameter, typically (but not necessarily) in
   \eqn{(0,1]}.  The coefficient \eqn{\lambda} of the integral of the
   squared second derivative in the fit (penalized log likelihood)
   criterion is a monotone function of \code{spar}, see the details
   below.}
 \item{cv}{ordinary (\code{TRUE}) or \dQuote{generalized} cross-validation
   (GCV) when \code{FALSE}.}
 \item{all.knots}{if \code{TRUE}, all distinct points in \code{x} are used as
   knots.  If \code{FALSE} (default), a subset of \code{x[]} is used,
   specifically \code{x[j]} where the \code{nknots} indices are evenly
   spaced in \code{1:n}, see also the next argument \code{nknots}.}
 \item{nknots}{integer giving the number of knots to use when
   \code{all.knots=FALSE}.  Per default, this is less than \eqn{n}, the
   number of unique \code{x} values for \eqn{n > 49}.}
 \item{df.offset}{allows the degrees of freedom to be increased by
   \code{df.offset} in the GCV criterion.}
 \item{penalty}{the coefficient of the penalty for degrees of freedom
   in the GCV criterion.}
 \item{control.spar}{optional list with named components controlling the
   root finding when the smoothing parameter \code{spar} is computed,
   i.e., missing or \code{NULL}, see below.

   \bold{Note} that this is partly \emph{experimental} and may change
   with general spar computation improvements!

   \describe{
     \item{low:}{lower bound for \code{spar}; defaults to -1.5 (used to
       implicitly default to 0 in \R versions earlier than 1.4).}
     \item{high:}{upper bound for \code{spar}; defaults to +1.5.}
     \item{tol:}{the absolute precision (\bold{tol}erance) used; defaults
     to 1e-4 (formerly 1e-3).}
     \item{eps:}{the relative precision used; defaults to 2e-8 (formerly
       0.00244).}
     \item{trace:}{logical indicating if iterations should be traced.}
     \item{maxit:}{integer giving the maximal number of iterations;
       defaults to 500.}
   }
   Note that \code{spar} is only searched for in the interval
   \eqn{[low, high]}.
 }
}
\details{
  The \code{x} vector should contain at least four distinct values.
  \emph{Distinct} here means \dQuote{distinct after rounding to 6 significant
  digits}, i.e., \code{x} will be transformed to
  \code{unique(sort(signif(x, 6)))}, and \code{y} and \code{w} are
  pooled accordingly.

  The computational \eqn{\lambda} used (as a function of
  \eqn{s=spar}{\code{spar}}) is
  \eqn{\lambda = r * 256^{3 s - 1}}{lambda = r * 256^(3*spar - 1)}
  where
  \eqn{r = tr(X' W X) / tr(\Sigma)},
  \eqn{\Sigma} is the matrix given by
  \eqn{\Sigma_{ij} = \int B_i''(t) B_j''(t) dt}{%
    Sigma[i,j] = Integral B''[i](t) B''[j](t) dt},
  \eqn{X} is given by \eqn{X_{ij} = B_j(x_i)}{X[i,j] = B[j](x[i])},
  \eqn{W} is the diagonal matrix of weights (scaled such that
  its trace is \eqn{n}, the original number of observations)
  and \eqn{B_k(.)}{B[k](.)} is the \eqn{k}-th B-spline.

  Note that with these definitions, \eqn{f_i = f(x_i)}, and the B-spline
  basis representation \eqn{f = X c} (i.e., \eqn{c} is
  the vector of spline coefficients), the penalized log likelihood is
  \eqn{L = (y - f)' W (y - f) + \lambda c' \Sigma c}, and hence
  \eqn{c} is the solution of the (ridge regression)
  \eqn{(X' W X + \lambda \Sigma) c = X' W y}.

  If \code{spar} is missing or \code{NULL}, the value of \code{df} is used to
  determine the degree of smoothing.  If both are missing, leave-one-out
  cross-validation (ordinary or \dQuote{generalized} as determined by
  \code{cv}) is used to determine \eqn{\lambda}.
  Note that from the above relation,
%%  lam      = r * 256^(3s - 1)
%%  log(lam) = log(r) + (3s - 1) * log(256)
%% (log(lam) - log(r)) / log(256)  = 3s - 1
%% s = [1 +  {log(lam) - log(r)} / {8 log(2)} ] / 3
%%   = 1/3 + {log(lam) - log(r)} / {24 log(2)}
%%   = 1/3 - log(r)/{24 log(2)} +  log(lam) / {24 log(2)}
%%   =               s0         + 0.0601 * log(lam)
  \code{spar} is \eqn{s = s0 + 0.0601 * \bold{\log}\lambda}{%
                   spar = s0 + 0.0601 * log(lambda)},
  which is intentionally \emph{different} from the S-plus implementation
  of \code{smooth.spline} (where \code{spar} is proportional to
  \eqn{\lambda}).  In \R's (\eqn{\log \lambda}) scale, it makes more
  sense to vary \code{spar} linearly.

  Note however that currently the results may become very unreliable
  for \code{spar} values smaller than about -1 or -2.  The same may
  happen for values larger than 2 or so. Don't think of setting
  \code{spar} or the controls \code{low} and \code{high} outside such a
  safe range, unless you know what you are doing!

  The \dQuote{generalized} cross-validation method will work correctly when
  there are duplicated points in \code{x}.  However, it is ambiguous what
  leave-one-out cross-validation means with duplicated points, and the
  internal code uses an approximation that involves leaving out groups
  of duplicated points.  \code{cv=TRUE} is best avoided in that case.
}
\note{
  The default \code{all.knots = FALSE} and \code{nknots = NULL} entails
  using only \eqn{O(n^{0.2})}
  knots instead of \eqn{n} for \eqn{n > 49}.  This cuts speed and memory
  requirements, but not drastically anymore since \R version 1.5.1 where
  it is only \eqn{O(n_k) + O(n)}{O(nk) + O(n)} where \eqn{n_k}{nk} is
  the number of knots.
  In this case where not all unique \code{x} values are
  used as knots, the result is not a smoothing spline in the strict
  sense, but very close unless a small smoothing parameter (or large
  \code{df}) is used.
}
\value{
  An object of class \code{"smooth.spline"} with components
  \item{x}{the \emph{distinct} \code{x} values in increasing order, see
  the \bold{Details} above.}
  \item{y}{the fitted values corresponding to \code{x}.}
  \item{w}{the weights used at the unique values of \code{x}.}
  \item{yin}{the y values used at the unique \code{y} values.}
  \item{lev}{leverages, the diagonal values of the smoother matrix.}
  \item{cv.crit}{cross-validation score, \dQuote{generalized} or true, depending
    on \code{cv}.}
  \item{pen.crit}{penalized criterion}
  \item{crit}{the criterion value minimized in the underlying
    \code{.Fortran} routine \file{sslvrg}.}
  \item{df}{equivalent degrees of freedom used.  Note that (currently)
    this value may become quite unprecise when the true \code{df} is
    between and 1 and 2.
  }
  \item{spar}{the value of \code{spar} computed or given.}
  \item{lambda}{the value of \eqn{\lambda} corresponding to \code{spar},
    see the details above.}
  \item{iparms}{named integer(3) vector where \code{..$ipars["iter"]}
    gives number of spar computing iterations used.}
  \item{fit}{list for use by \code{\link{predict.smooth.spline}}, with
    components
    \describe{
      \item{knot:}{the knot sequence (including the repeated boundary
	knots).}
      \item{nk:}{number of coefficients or number of \dQuote{proper}
	knots plus 2.}
      \item{coef:}{coefficients for the spline basis used.}
      \item{min, range:}{numbers giving the corresponding quantities of
	\code{x}.}
    }
  }
  \item{call}{the matched call.}
}
\references{
  Chambers, J. M. and Hastie, T. J. (1992)
  \emph{Statistical Models in S}, Wadsworth \& Brooks/Cole.

  Green, P. J. and Silverman, B. W. (1994)
  \emph{Nonparametric Regression and Generalized Linear Models:
    A Roughness Penalty Approach.} Chapman and Hall.

  Hastie, T. J. and Tibshirani, R. J. (1990)
  \emph{Generalized Additive Models.}  Chapman and Hall.
}
\author{
  \R implementation by B. D. Ripley and Martin Maechler
  (\code{spar/lambda}, etc).

  This function is based on code in the \code{GAMFIT} Fortran program by
  T. Hastie and R. Tibshirani (\url{http://lib.stat.cmu.edu/general/}),
  which makes use of spline code by Finbarr O'Sullivan.  Its design
  parallels the \code{smooth.spline} function of Chambers \& Hastie (1992).
}
\seealso{
  \code{\link{predict.smooth.spline}} for evaluating the spline
  and its derivatives.
}
\examples{
data(cars)
attach(cars)
plot(speed, dist, main = "data(cars)  &  smoothing splines")
cars.spl <- smooth.spline(speed, dist)
(cars.spl)
## This example has duplicate points, so avoid cv=TRUE
\dontshow{
  all(cars.spl $ w == table(speed)) # TRUE (weights = multiplicities)
  str(cars.spl, digits=5, vec.len=6)
  cars.spl$fit
}
lines(cars.spl, col = "blue")
lines(smooth.spline(speed, dist, df=10), lty=2, col = "red")
legend(5,120,c(paste("default [C.V.] => df =",round(cars.spl$df,1)),
               "s( * , df = 10)"), col = c("blue","red"), lty = 1:2,
       bg='bisque')
detach()

##-- artificial example
y18 <- c(1:3,5,4,7:3,2*(2:5),rep(10,4))
xx  <- seq(1,length(y18), len=201)
(s2  <- smooth.spline(y18)) # GCV
(s02 <- smooth.spline(y18, spar = 0.2))
plot(y18, main=deparse(s2$call), col.main=2)
lines(s2, col = "gray"); lines(predict(s2, xx), col = 2)
lines(predict(s02, xx), col = 3); mtext(deparse(s02$call), col = 3)

## The following shows the problematic behavior of 'spar' searching:
(s2  <- smooth.spline(y18,          con=list(trace=TRUE,tol=1e-6, low= -1.5)))
(s2m <- smooth.spline(y18, cv=TRUE, con=list(trace=TRUE,tol=1e-6, low= -1.5)))
## both above do quite similarly (Df = 8.5 +- 0.2)

\dontshow{
% not for the user
  s2. <- smooth.spline(y18, cv=TRUE,con=list(trace=TRUE, tol=1e-6,low= -3,maxit=20))
  s2. ## Intel-Linux: Df ~= (even! > ) 18 : interpolating -- much smaller PRESS
     ## {others, e.g., may end quite differently!}
  lines(predict(s2., xx), col = 4)
  mtext(deparse(s2.$call,200), side= 1, line= -1, cex= 0.8, col= 4)

  sdf8 <- smooth.spline(y18, df = 8, con=list(trace=TRUE))
  sdf8 ; sdf8$df - 8

  try(smooth.spline(y18, spar = 50)) #>> error : spar 'way too large'
}
}
\keyword{smooth}
\eof
\name{smoothEnds}
\title{End Points Smoothing (for Running Medians)}
\alias{smoothEnds}
\usage{
smoothEnds(y, k = 3)
}
\description{
  Smooth end points of a vector \code{y} using subsequently smaller
  medians and Tukey's end point rule at the very end. (of odd span),
}
\arguments{
  \item{y}{dependent variable to be smoothed (vector).}
  \item{k}{width of largest median window; must be odd.}
}
\value{vector of smoothed values, the same length as \code{y}.
}
\details{
  \code{smoothEnds} is used to only do the \dQuote{end point smoothing},
  i.e., change at most the observations closer to the beginning/end
  than half the window \code{k}.  The first and last value are computed using
  \dQuote{\emph{Tukey's end point rule}}, i.e.,
  \code{sm[1] = median(y[1], sm[2], 3*sm[2] - 2*sm[3])}.
}
\references{
  John W. Tukey (1977)
  \emph{Exploratory Data Analysis}, Addison.

  Velleman, P.F., and Hoaglin, D.C. (1981)
  \emph{ABC of EDA (Applications, Basics, and Computing of Exploratory
    Data Analysis)}; Duxbury.
}
\author{Martin Maechler}
\seealso{
  \code{\link{runmed}(*, end.rule = "median")} which calls
  \code{smoothEnds()}.
}
\examples{
y <- ys <- (-20:20)^2
y [c(1,10,21,41)] <-  c(100, 30, 400, 470)
s7k <- runmed(y,7, end = "keep")
s7. <- runmed(y,7, end = "const")
s7m <- runmed(y,7)
col3 <- c("midnightblue","blue","steelblue")
plot(y, main = "Running Medians -- runmed(*, k=7, end.rule = X)")
lines(ys, col = "light gray")
matlines(cbind(s7k,s7.,s7m), lwd= 1.5, lty = 1, col = col3)
legend(1,470, paste("end.rule",c("keep","constant","median"),sep=" = "),
       col = col3, lwd = 1.5, lty = 1)

stopifnot(identical(s7m, smoothEnds(s7k, 7)))
}
\keyword{smooth}
\keyword{robust}
\eof
\name{supsmu}
\alias{supsmu}
\title{Friedman's SuperSmoother}
\description{
  Smooth the (x, y) values by Friedman's \dQuote{super smoother}.
}
\usage{
supsmu(x, y, wt, span = "cv", periodic = FALSE, bass = 0)
}
\arguments{
  \item{x}{x values for smoothing}
  \item{y}{y values for smoothing}
  \item{wt}{case weights, by default all equal}
  \item{span}{the fraction of the observations in the span of the running
    lines smoother, or \code{"cv"} to choose this by leave-one-out
    cross-validation.}
  \item{periodic}{if \code{TRUE}, the x values are assumed to be in
    \code{[0, 1]} and of period 1.}
  \item{bass}{controls the smoothness of the fitted curve. Values of up
    to 10 indicate increasing smoothness.}
}
\details{
  \code{supsmu} is a running lines smoother which chooses between three
  spans for the lines. The running lines smoothers are symmetric, with
  \code{k/2} data points each side of the predicted point, and values of
  \code{k} as \code{0.5 * n}, \code{0.2 * n} and \code{0.05 * n}, where
  \code{n} is the number of data points.  If \code{span} is specified,
  a single smoother with span \code{span * n} is used.

  The best of the three smoothers is chosen by cross-validation for each
  prediction. The best spans are then smoothed by a running lines
  smoother and the final prediction chosen by linear interpolation. 
 
  The FORTRAN code says: \dQuote{For small samples (\code{n < 40}) or if
    there are substantial serial correlations between observations close
    in x-value, then a prespecified fixed span smoother (\code{span >
      0}) should be used.  Reasonable span values are 0.2 to 0.4.}
}
\value{
  A list with components
  \item{x}{the input values in increasing order with duplicates removed.}
  \item{y}{the corresponding y values on the fitted curve.}
}
\references{
  Friedman, J. H. (1984)
  SMART User's Guide.
  Laboratory for Computational Statistics, Stanford University Technical
  Report No. 1.

  Friedman, J. H. (1984)
  A variable span scatterplot smoother.
  Laboratory for Computational Statistics, Stanford University Technical
  Report No. 5.
}
\seealso{\code{\link{ppr}}}

\examples{
data(cars)
with(cars, {
    plot(speed, dist)
    lines(supsmu(speed, dist))
    lines(supsmu(speed, dist, bass = 7), lty = 2)
    })
}
\keyword{smooth}
\eof
