Spline fitting in R - how to force passing two data points?

Question

I am using "smooth.spline" in R. Here is a snippet from the documentation:

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/smooth.spline.html

smooth.spline {stats} R Documentation

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, keep.data = TRUE, df.offset = 0, penalty = 1, control.spar = list(), tol = 1e-6 * IQR(x))

---

My question is that:

I have two vectors of data x and y, where the lower bound for x is -100 and the upper bound for x is +100.

And I knew that for y=f(x):

f(-100)=-1

f(+100)=+1

That's to say, I would like to have the lower boundary and upper boundary points to be forced passing thru by the cubic spline procedure.

Because these two points are accurate and precise.

How to do that?

Could anybody please help me?

Thanks a lot!

Answer 1

Rather than use smooth.spline() in the stats package, there is a function cobs() in the cobs package that allows you to do exactly the sort of thing you want. COBS stands for Constrained B-splines. Possible constraints include going through specific points, setting derivatives to specified values, monotonicity (increasing or decreasing), concavity, convexity, periodicity, etc.

In your case, use

cobs(x, y, pointwise=rbind(c(0,-100,-1),c(0,100,1)))

Answer 2

I cannot think of any way to do it using smooth.spline. If you were to use a spline basis such as bs from the splines package, then you could possibly do this using quadradic programming to constrain the endpoints, but it could be complicated figuring out the constraints.

Here is an approach that uses xsplines (different but similar to other types of splines) and the optim function to find the values to use (nls could be used as well). I chose 3 internal equally spaced control points and a shape of 1, but you could play with these to compare the fit:

x <- seq( -100, 100, length=101 )
y <- sin( x/200*pi ) + rnorm(101, 0, 0.15)

myfun <- function(par) {
    yh <- c(-1, par, 1)
    xh <- c(-100, -50, 0, 50, 100)
    sp <- xspline( xh,yh, shape=1, draw=FALSE)
    yhat <- approx( sp, xout=x )$y
    sum( (y-yhat)^2 )
}

out <- optim( c(-.5, 0, .5), myfun )

plot(x,y)
xspline( c(-100, -50, 0, 50, 100), c(-1, out$par, 1), shape=1,
        border='blue' )