Interpolating Surface in R

Question

I'm trying to interpolate a polynomial of the form

z = Ax^2 + By^2 + Cxy + Dx + Ey + F

In R, where the capital letters are constant coefficients.

For the following data

My horizontal axis is: KSo<-c(0.90,0.95,1.00,1.05,1.10)

My vertical axis is: T<-c(1/12,3/12,6/12,1,2,5)

And the data mapped by KSo X T are:

14.2 13.0 12.0 13.1 14.5

14.0 13.0 12.0 13.1 14.2

14.1 13.3 12.5 13.4 14.3

14.7 14.0 13.5 14.0 14.8

15.0 14.4 14.0 14.5 15.1

14.8 14.6 14.4 14.7 15.0

In other words, the observed datum for (1.00,6/12) is 12.5

How would I interpolate, for example, the predicted data for (0.98,11/12)?

Edit: I found a nice package, akima, with the bicubic function, that uses splines. I'd still like to see what people suggest

Answer 1

Here is the two possible suggested solutions.
The original data:

KSo<-c(0.90,0.95,1.00,1.05,1.10)
T<-c(1/12,3/12,6/12,1,2,5)
mapping<-c(14.2, 13.0, 12.0, 13.1, 14.5,
           14.0, 13.0, 12.0, 13.1, 14.2,
           14.1, 13.3, 12.5, 13.4, 14.3,
           14.7, 14.0, 13.5, 14.0, 14.8,
           15.0, 14.4, 14.0, 14.5, 15.1,
           14.8, 14.6, 14.4, 14.7, 15.0)
mapped<-matrix(mapping, ncol=5, byrow=TRUE)

Here is the linear interpolation solution:

#predict
x<-0.98
y<-11/12

#Perform 2D interpolation
#find index along x and y axis
ki<-as.integer(cut(x, KSo, right=FALSE))
Ti<-as.integer(cut(y, T, right=FALSE))

#dx = (x-x1)/(x2-x1)  where x is the point to interpolate to.
# dx will vary from 0 to <1, (or worded differently the % distance between the 2 grid points)
dx<-(x-KSo[ki])/(KSo[ki+1]-KSo[ki])
dy<-(y-T[Ti])/(T[Ti+1]-T[Ti])

#find values and weighed sums of neighboring points
#    equations as per Wikipedia
f00<-mapped[Ti, ki]*(1-dx)*(1-dy)  
#    f(0,0) -weight each corner contributes to the final results
f01<-mapped[Ti+1, ki]*(1-dx)*dy
f10<-mapped[Ti, ki+1]*dx*(1-dy)
f11<-mapped[Ti+1, ki+1]*dx*dy
sum(f00, f10, f01, f11)

Same analysis as above but with R's functions

ki<-as.integer(cut(x, KSo, right=FALSE))
Ti<-as.integer(cut(y, T, right=FALSE))
ilower<-approx(T, mapped[,ki], y)$y
iupper<-approx(T, mapped[,(ki+1)], y)$y
approx(KSo[ki:(ki+1)], c(ilower, iupper), x)

Below is the Regression Model using the entire data set to fit. Since the data is not perfect fit, the estimated value at a grid point is different than the original specified value.

#establish data frame with raw data
df<-data.frame(expand.grid(KSo, T), mapping)
names(df)<-c("KSo", "T", "z")

#Perform regression
model<-lm(z~I(KSo^2)+I(T^2)+KSo*T, data=df)

#Predict results at desired coordinates
dfnew<-data.frame(KSo=c(1, 0.98), T=c(0.5, 0.9166667))
predict(model, dfnew)

Hope this helps.