Ch4.2: Piecewise Interpolation

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

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

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

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

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

• 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

• 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

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

x = [-2, -1, 0, 1]
y = [-1, -2, -1, 2]
p = interp1d(x, y, kind = 'cubic')

n = len(x)
xnew = np.linspace(x[0], x[n-1], num=500, endpoint=True)
plt.plot(x, y, 'o', xnew, p(xnew), '-,r')
plt.legend(['Original Data','Cubic Spline p(x)'], loc = 'best')
plt.show()

• Wage and other data for 3000 male workers in Mid-Atlantic region; see weblink.

#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])

• Extra slides for background programs

polyinterp <- function(x,y) {
  if(length(x) != length(y))
   stop("Length of x & y vectors must be same")
  n <- length(x) - 1
  vandermonde <- rep(1,length(x))
  for(i in 1:n) {
    xi <- x^i
    vandermonde <- cbind(vandermonde,xi) }
  beta <- solve(vandermonde,y)
  names(beta) <- NULL
  return(beta) }

• Nested polynomial evaluation

\[ \begin{align*} f(x) & = 1 - 2x + 3x^2 +4x^3 +5x^4 \\ & = 1 + x(-2 + x(3 + x(4 + 5x)) \end{align*} \]

horner <- function(x, coefs) {
  y <- rep(0, length(x))
  for(i in length(coefs):1)
    y <- coefs[i] + x*y
  return (y)  }

HornerPlot <- function(x,y) {
  n = length(x)
  N <- 100; h <- (x[n] - x[1])/N
  t <- rep(0, N); f <- rep(0, N) 
  t[1] <- x[1]; f[1] <- y[1]; 
  p <- polyinterp(x, y)

  for(i in 1:N) { 
     t[i+1] <- t[i] + h 
     f[i+1] <- horner(t[i+1],p) }

  plot(t,f, 
    main = "Horner Interpolating Polynomial",
    xlab="x",ylab="y",type="l",lwd=2,col="red", 
    xlim = c(x[1],x[n]), ylim = c(min(f),max(f)) )
  lines(x,y,type = "p",col="blue",lwd=3)}