Ch4.2: Piecewise Interpolation

Higher Degree Polynomial Interpolation

  • Suppose we had a data set of 100 data points.
  • A polynomial interpolation to this data would have degree 99.
  • Can have as much as 99 oscillations, typically of large magnitude.
  • These are found by taking the derivative (deg 99 polynomial) and setting equal to zero, and then solving for critical points.

title

Piecewise Interpolation

  • A piecewise interpolation uses lower order polynomials between successive portions of data.
  • Common choices are linear and cubic.
  • Typically referred to as splines.

title

Quadratic Spline Figure

  • Note the piecewise listing of quadratic formulas

title

Draftsman's Spline

title

  • Woodworking
  • Ship hull design
  • Fit flexible strip of wood through pegs

Clamped Spline

title

  • Spline is clamped at endpoints
  • Specific derivative (shape) condition at endpoints

Natural Spline

title

  • Spline is not clamped at endpoints (free)
  • Linear shape at endpoints
  • Second derivative equal to zero at endpoints.

Linear Interpolation with R

  • Program from our book.
  • Fit equation of line to two data points \( (x_1,y_1), \, (x_2,y_2) \)
  • For \( y = mx + b \), need \( m \) and \( b \)
linterp <- function (x1, y1, x2, y2) {
  m <- (y2-y1)/(x2-x1)
  b <- y2 - m*x2
  ## Convert into a form suitable for horner
  return (c(b, m))  }

Linear Spline with R

  • Program from our book.
  • Returns coefficients for each linear spline. 
pwiselinterp <- function (x, y) {
  n <- length(x) - 1
  y <- y[order(x)]
  x <- x[order(x)]
  mvec <- bvec <- c()
  for(i in 1:n) {
    p <- linterp(x[i], y[i], x[i + 1], y[i + 1])
    mvec <- c(mvec,p[2])
    bvec <- c(bvec,p[1]) }
  return(list(m = mvec , b = bvec )) }

Example 1: Linear Spline

x <- c(-2, -1, 0, 1)
y <- c(-1, -2, -1, 2)
pwiselinterp(x, y)

$m
[1] -1  1  3

$b
[1] -3 -1 -1

Data Plot Linear Function 

DataPlotLinear <- function(x,y) {
  n = length(x)
  xplotmin <- x[1]-1
  xplotmax <- x[n]+1
  yplotmin <- min(c(min(y),0))-1
  yplotmax <- max(y)+1

  plot(x,y, 
    main = "Data Plot",
    xlab="x",ylab="y",type="p",lwd=3,col="blue", 
    xlim = c(xplotmin,xplotmax), 
    ylim = c(yplotmin,yplotmax) ) 
  lines(x,y,type = "l",col="red",lwd=2)}

Example 1: Linear Spline

x <- c(-2, -1, 0, 1)
y <- c(-1, -2, -1, 2)
pwiselinterp(x, y)

$m
[1] -1  1  3

$b
[1] -3 -1 -1

DataPlotLinear(x,y)

plot of chunk unnamed-chunk-6

Example 1: R Function

  • R has a piecewise interpolation function
x <- c(-2, -1, 0, 1)
y <- c(-1, -2, -1, 2)
f <- approxfun(x, y)
f(0.8)

[1] 1.4

Example 1: R Function & Plot

  • Using R's piecewise interpolation function
x <- c(-2, -1, 0, 1)
y <- c(-1, -2, -1, 2)
f <- approxfun(x, y)
f(0.8)

[1] 1.4

DataPlotLinear(x,y)

plot of chunk unnamed-chunk-9

Example 2: R Cubic Spline Function

  • R has a cubic spline function
  • Returns interpolated y-values at specified x-values
x <- c(-2, -1, 0, 1)
y <- c(-1, -2, -1, 2)
spline.ex <- splinefun(x, y, method = "natural")
spline.ex(x)

[1] -1 -2 -1  2

spline.ex(0.8)

[1] 1.3232

Cubic Spline Plot using R

library(pracma)

CubicSplinePlot <- function(x,y) {
  n = length(x);  N <- 100 
  t <- linspace(x[1], x[n], N) 
  p <- cubicspline(x, y, t)
  bottom = min(0,min(p)); top = max(0,max(p))
plot(t,p, 
   main = "Cubic Spline",
   xlab="x",ylab="y",type="l",lwd=2,col="red", 
   xlim = c(x[1],x[n]), ylim = c(bottom, top) )
points(x,y,type = "p",col="blue",lwd=3)
legend("topleft",
   legend = c("Data", "Cubic Spline"),
   col=c("blue","red"),  
   lty=c(1,1) )}

Example 2: Cubic Spline FUnction & Plot

x <- c(-2, -1, 0, 1)
y <- c(-1, -2, -1, 2)
spline.ex(0.8)

[1] 1.3232

CubicSplinePlot(x,y)

plot of chunk unnamed-chunk-14

Example 3: Wage Data

  • Wage and other data for 3000 male workers in Mid-Atlantic region; see weblink.
  • Uses ISLR package. 
#loading the Splines Packages
require(splines)
#ISLR contains the Dataset
require(ISLR)
attach(Wage) #attaching Wage dataset
#?Wage #for more details on the dataset
agelims<-range(age)
#Generating Test Data
age.grid<-seq(from=agelims[1], to = agelims[2])

Example 3: Wage Data

  • Additional commands
#3 cutpoints at ages 25 ,50 ,60
fit<-lm(wage ~ bs(age,knots = c(25,40,60)),data = Wage )
#summary(fit)
#Plotting the Regression Line to the scatterplot   
plot(age,wage,col="grey",xlab="Age",ylab="Wages")
points(age.grid,predict(fit,newdata = list(age=age.grid)),col="darkgreen",lwd=2,type="l")
#adding cutpoints
abline(v=c(25,40,60),lty=2,col="darkgreen")

Example 3: Wage Data

  • Output plot

plot of chunk unnamed-chunk-17