Ch4.2: Piecewise Interpolation

What do you call a fake noodle?

An impasta.

I like telling Dad jokes.

Sometimes he laughs!

• 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 \).

\[ \small{ \begin{aligned} m & = \frac{y_2 - y_1}{x_2-x_1} \\ b & = y_2 - m x_2 \end{aligned} } \]

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))  }

• 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 )) }

• Graphs the linear spline with data

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 = "Linear Splines",
    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)}

• Instead of trying to fit a piecewise linear model through a set of data points, we can try to fit a piecewise cubic polynomial model through the points.
  
• For a smooth curve (no jaggedness), we require that the resulting graph is differentiable everywhere.

• Suppose that we want to fit cubic splines through the data points (1,2), (2,1), (3,2).
  
• Find two cubic polynomials (splines) that fit their respective two data points along with additional smoothness conditions.

\[ \begin{aligned} p_1(x) &= a_1 x^3 + b_1 x^2 + c_1 x + d_1 , \,\, 1 \leq x \leq 2 \\ p_2(x) &= a_2 x^3 + b_2 x^2 + c_2 x + d_2 , \,\, 2 \leq x \leq 3 \end{aligned} \]

• Here, \( p_1(x) \) passes through (1,2) and (2,1).
• Also, \( p_2(x) \) passes through (2,1) and (3,2).
• To join smoothly at the shared point (2,1), they need to have the same slope and concavity there.

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 Splines",
   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)}

• This application turned up on an internet search.
  
• 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])

• 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")

For our next class period, you will need a computer that runs RStudio and a recent version of R.

• This could be one of the classroom computers, or a laptop.
  
• We will also need to install and use the pracma package.

For the project, we collect data and use R to fit a cubic spline.