curfit: Curve fitting with splines

Description Usage Arguments Details Value Author(s) References See Also Examples

Description

General curve fitting with splines. Wrapper for Fortran function CURFIT.

Usage

1
2
3
4
5
6
curfit(x, ...)
percur(x, ...)
## Default S3 method:
curfit(x, y = NULL, w = NULL, s=NULL,
       knots = NULL, n = NULL, from = min(x), to = max(x),
       k = 3, periodic = FALSE, ...)

Arguments

x

A data.frame, matrix, or numeric vector. See details.

y

Optional numeric vector. See details.

w

Optional vector of weights

s

a nonnegative number of NULL.

knots

Optional vector of knots including end knots. See details.

n

The number of knots. If both knots and n are provided, n must equal length(knots).

from

Lower bound from which to fit spline.

to

Upper bound to which to fit spline.

k

Degree of the spline. Valid options for k are 1 to 5, inclusively.

periodic

logical; if TRUE then s(a) == s(b)

...

Additional arguments used only in update.curfit. Otherwise, ignored.

Details

curfit determines a spline approximation s(x) of degree k (order = k+1). See Dierckx (1993, ch. 1 and 3) for definition of symbols.

curfit uses two alternative methods: least squares an smoothing spline.

As with smooth.spline, the x vector should contain at least four distinct values. Distinct here means “distinct after rounding to 6 significant digits”, i.e., x will be transformed to unique(signif(x, 6)), and y and w are pooled accordingly.

NOTE: curfit.f calls fpchec.f, which checks to ensure that the number of knots n is at least twice the order of the spline, i.e., 2*(k+1), where k is the degree of the spline. It also checks that the number of distinct data points exceeds the degree k of the spline.

For the default method, arguments x, y and w are supplied to xyw.coords to determine abscissa and ordinate values. If any non-distinct x values are found, they are combined, replacing the corresponding values of y and w by the mean and sum, respectively.

When supplying knots, the end knots must be included. Thus, if periodic = FALSE, the first k+1 values should be min(x) and the last n+k-1 value should be max(x). See below for examples.

If neither knots nor n are provided, the algorithm places one knot at each distinct x value. If n is provided but not knots, n-2*(k+1) interior knots are evenly spaced between from and to. An error message is given if both knots and n are provided and n does not equal length(knots).

percur is equivalent to calling curfit(..., periodic=TRUE).

Value

An object of class dierckx with the following components:

iopt

method used coded as follows:

-1

least squares using all knots

0

smoothing spline with a subset of knots

1

smoothing spline update from a previous call

m

length of 'x'

x

abscissa values

y

ordinate values

w

input weights

from

input from value

to

input to value

k

degree of the spline

s

input smoothing parameter

nest

Estimated number of knots

n

Actual number of knots

knots

Knot locations. Use knots.dierckx to extract.

coef

B-spline coefficients. Use coef.dierckx to extract.

fp

sum of squares residuals. Use deviance.dierckx to extract.

wrk

Work space. Do NOT modify before call to update.dierckx.

lwrk

Length of wrk. Do NOT modify before call to update.dierckx.

iwrk

Integer work space. Do NOT modify before call to update.dierckx.

ier

Error code. Should always be zero.

message

brief character string description of the fit

g

Number of interior knots

method

method coded 'ls' for least squares (if 's' is not provided), or 'ss' for smoothing spline (if 's' is provided).

periodic

input 'periodic' parameter

routine

Always 'curfit.default'

xlab

The x-label determined from deparse(substitute(x)).

ylab

The y-label determined from deparse(substitute(y)).

Author(s)

Sundar Dorai-Raj and Spencer Graves

References

Dierckx, P. (1993) Curve and Surface Fitting with Splines, Oxford Science Publications.

See Also

concon, spline, smooth.spline

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
##
## titanium example
##
data(titanium)
# The following kills R:  
# with(titanium, curfit(x, y))


##
## made up example
##
x <- 0:24
y <- c(1.0,1.0,1.4,1.1,1.0,1.0,4.0,9.0,13.0,
       13.4,12.8,13.1,13.0,14.0,13.0,13.5,
       10.0,2.0,3.0,2.5,2.5,2.5,3.0,4.0,3.5)

#fitLS0 <- curfit(x, y)
#fitSS0 <- curfit(x, y, s=0)

ks <- c(3, 5, 2)
kk <- length(ks) 
z <- vector("list", kk) 
names(z) <- ks 
for(i in 1:kk) {
  k <- ks[i] 
  z1 <- curfit(x, y, s = 1000, k = k)
  z2 <- update(z1, s = 60)
  z3 <- update(z2, s = 10)
  z4 <- update(z3, s = 30)
  z5 <- curfit(x, y, s = 30, k = k)
  z6 <- update(z5, s = 0)
  knots <- c(rep(0, k + 1), seq(3, 21, 3), rep(24, k + 1))
  z7 <- curfit(x, y, s = 30, knots = knots, k = k)
  z[[i]] <- list(z1, z2, z3, z4, z5, z6, z7)
}

p <- unlist(z, recursive = FALSE)
n <- sapply(lapply(p, knots), length)
s <- sapply(p, "[[", "s")
i <- sapply(p, "[[", "iopt")
m <- ifelse(i == -1, "ls", ifelse(i == 0, "ss", "ss1"))
k <- sprintf("k = %d", sapply(p, "[[", "k"))
g <- sprintf("%s(s=%d)", m, s, i)
sp <- data.frame(x = rep(x, times = length(p)),
     y = rep(y, times = length(p)), z = unlist(lapply(p, fitted)),
     k = factor(rep(k, each = length(x))), g = rep(g, each = length(x)))

library(lattice)
xyplot(z ~ x | k, data = sp, groups = g,
   panel = function(x, y, subscripts, groups, obs, ...) {
     panel.superpose(x, y, subscripts, groups, lwd = 3, type = "l", ...)
     x <- unique(x)
     y <- unique(obs)
     panel.xyplot(x, obs, pch = 16, cex = 1.2, col = "darkblue")
   },
   auto.key = list(space = "right", points = FALSE, lines = TRUE),
   obs = sp$y)

## periodic spline
set.seed(42)
n <- 100
r <- 1:n
x <- 0.01 * (r - 1)
e <- rnorm(n, 0, 0.1)
w <- rep(1/sd(e), n + 1)
y <- cos(2 * pi * x) + 0.25 * sin(8 * pi * x) + e
x <- c(x, 1)
y <- c(y, y[1])
kn <- seq(0.01, 0.99, length = 12)
f1 <- percur(x, y, w = w, s = 90, k = 5)

library(lattice)
top <- xyplot(y ~ x,
              panel = function(x, y, ...) {
                panel.abline(v = knots(f1), lty = 2, lwd = 3, col = "gray")
                panel.xyplot(x, y, pch = 16, col = "#800000", cex = 1.2)
                panel.xyplot(x, fitted(f1), type = "l", lwd = 3, col = "#000080")
              },
              par.settings = list(layout.widths = list(left.padding = 0, right.padding = 0)),
              scales = list(cex = 1.2),
              xlab = "", ylab = "")
newx <- seq(-2, 2, 0.01)
newy <- predict(f1, newx)
bot <- xyplot(newy ~ newx, type = "l",
              panel = function(...) {
                panel.abline(v = -2:2, lty = 2, col = "salmon", lwd = 3)
                panel.xyplot(...)
              },
              col = "#000080", lwd = 3,
              par.settings = list(layout.widths = list(left.padding = 0, right.padding = 0)),
              scales = list(cex = 1.2),
              xlab = "", ylab = "")
print(top, c(0, 0.2, 1, 1))
print(bot, c(0.008, 0, 0.992, 0.25), newpage = FALSE)

## example borrowed from ?smooth.spline
plot(cars$speed, cars$dist,
     main = "data(cars)  &  smoothing splines",
     xlab = "SPEED", ylab = "DISTANCE",
     cex.lab = 1.2, cex.axis = 1.2,
     cex.main = 2, cex = 1.5, col = "blue")
## This example has duplicate points, so avoid cv=TRUE
cars.spl.0 <- smooth.spline(cars$speed, cars$dist)
cars.spl.1 <- smooth.spline(cars$speed, cars$dist, df = 10)
cars.spl.2 <- curfit(cars$speed, cars$dist, s = 5e3)
newx <- seq(min(cars$speed), max(cars$speed), len = 200)
lines(predict(cars.spl.0, newx), col = "blue", lwd = 3, lty = 2)
lines(predict(cars.spl.1, newx), lty="dashed", col = "red", lwd = 3)
lines(newx, predict(cars.spl.2, newx), lty="dotted", lwd = 3)
legend(5, 120, c(paste("smooth.spline( * , df = ", round(cars.spl.0$df, 1), ")", sep = ""),
                 "smooth.spline( * , df = 10)", "curfit( * , s = 5e3)"),
       col = c("blue", "red", "black"),
       lty = c("solid", "dashed", "dotted"), lwd = 3,
       bg = 'bisque', cex = 1.5)

DierckxSpline documentation built on May 2, 2019, 6:30 p.m.