[This article was first published on plausibel, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
As a first step to produce some useable code for spline interpolation/approximation in R, I set out to first do polynomial interpolation to see how I get along. It’s not that there is no spline interpolation software for R, but I find it a bit limited. splinefun, for example, can do only 1-dimensional interpolation. interp{akima} can do bicubic splines – i.e. it can fit a spline surface to a set of points over gridded data as for example (x,y,f(x,y)). However, x and y need to have the same length. There’s no apparent reason for this, and it might be inefficient to have to do this at times, so I’ll be looking at changing that.
As a first step, I’ll implement a polynomial approximation algorithm. The basic principle is that of a convex hull of points in a given space. By repeatedly taking convex combinations of points one builds up the approximating curve. The theoretical background for this can be looked up here. The algorithm I implement here is a variant of the DeBoor triangular calculation. I am not claiming this to be the most efficient calculation, this is just a first try. Way too many loops.
So: I sample different number of points on a half circle and construct an interpolating polynomial. you can see the result in the graph below. the red line is the polynomial approximation. the code is added below the graph.
here is the code:
# test polynomial interpolation code by interpolating points on a semi-circle
# using a different number of points each time.
# contents:
# 1) make.circle(degree) is a function that generates
# n=degree+1 data points (x,y) on a semicircle. this is the "data" we
# want to interpolate.
# 2) pinterp() makes polynomial interpolation on "data" of degree "degree"
# using degree+1 parameters
# 3) run() is a wrapper that runs the interpolation
# note: polynomial interpolation uses n=degree+1 points and T=degree+1 parameters.
#
# output:
# plot of the interpolation overlaid by the used data points
rm(list=ls(all=TRUE))
setwd("~/dropbox/Rstuff/splines")
require(plotrix)
# this package contains a function to draw a circle.
# you could generate any other data as well.
# task: run this approximation for 3,7,11 and 15 sample points.
degrees <- c(2,6,10,14)
##################
# define functions
##################
make.circle <- function(degree){
# makes a vector of degree+1 points on a semi-circle
plot(-2:2,-2:2,type="n",xlab="",ylab="",main="make.circle output. disregard.")
vals <- draw.circle(0,0,radius=1,nv=2*(degree)+1)
data <- array(unlist(vals),dim=c(2*degree+1,2))
data <- data[data[ ,2]>=0, ] # take only upper half of circle§
return(data)
}
pinterp <- function(degree,params,data,t){
# polynomial interpolation on data at point t
# has a list "x" whose entries are each a
# matrix for points at t at approximation of degree
# deg. x gets shorter with increasing degree. at final entry
# degree+1, x is (1,2)
x <- vector("list",degree+1)
x[[1]] <- data # matrix 1 holds data for degree zero
# each row of x[[j]]
# represents (x,y) coords of a point
n <- degree+1 # number of points we're interpolating on
for (ideg in 2:n){ # recursively build interpolations of increasing degree
nn <- n-ideg+1
x[[ideg]] <- array(data=NA,dim=c(nn,2))
for (i in 1:nn){
x[[ideg]][i, ] <-
(params[i+ideg-1] - t)/(params[i+ideg-1] - params[i])*x[[ideg-1]][i, ] +
(t - params[i])/(params[i+ideg-1] - params[i])*x[[ideg-1]][i+1, ]
}
}
return(x[[degree+1]]) # return only final degree values,i.e. result
}
run <- function(degrees,t){
nd <- length(degrees)
rval <- vector("list",nd)
origdata <- vector("list",nd)
for (j in 1:nd){
degree <- degrees[j]
params <- seq(0,5,length=degree+1)
nt <- length(t)
rval[[j]] <- array(data=NA,dim=c(nt,2))
data <- make.circle(degree)
origdata[[j]] <- data
for (it in 1:nt) rval[[j]][it, ] <- pinterp(degree,params,data,t[it])
}
return(list(interpol=rval,orig=origdata))
}
###################
# produce output
###################
outrun <- run(degrees,t=seq(1,5,length=20))
rval <- outrun$interpol
orig <- outrun$orig
png(file="graphs/ex1.png")
par(mfrow=c(2,2))
for (j in 1:length(degrees)){
plot(x=orig[[j]][ ,1],y=orig[[j]][ ,2],type="o",
main=paste("degree: ",degrees[j],", points: ",
degrees[j]+1),xlab="x",ylab="y",ylim=c(0,1.5))
lines(x=rval[[j]][ ,1],y=rval[[j]][ ,2],col="red")
}
dev.off()
par(mfrow=c(1,1))flo.
To leave a comment for the author, please follow the link and comment on their blog: plausibel.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.