Polynomial interpolation

Generation of the raw data

n <- 15
set.seed(7654321)
x <- runif(n) # generate n random numbers in [0,1]
y <- runif(n, min = -1, 1) # generate response in [-1,1]
plot(x,y)

Polynomial interpolation of the raw data

interp <- lm(y ~ poly(x,n-1)) # polynomial interpolation
plot(x,y)
sx <- sort(x)
lines(sx, predict(interp, data.frame(x = sx)))

Prediction of a new data point

Generate a new data point

newx <- runif(1) 
newy <- runif(1) # add a new data point

Prediction

plot(c(0,1),c(-10,1), type = 'n', xlab = 'x', ylab = 'y')
points(x,y)
points(newx,newy, pch = 19)
lines(sort(c(x,newx)), predict(interp, data.frame(x = sort(c(x,newx)))))
points(sort(c(x,newx)), predict(interp, data.frame(x = sort(c(x,newx)))), pch = 3)

The difference between the data point and the prediction is huge.

print(predict(interp, data.frame(x = newx)) - newy)

Graph of the fitted polynomial

plot(c(0,1),c(-100,100), type = 'n', xlab = 'x', ylab = 'y')
points(x,y)
lines(seq(0,1, by = 0.005), predict(interp, data.frame(x = seq(0,1, by = 0.005))))

Linear regression on the same dataset

linreg <- lm(y ~ x) # linear regression
plot(x,y)
points(newx,newy, pch = 19)
sx <- sort(x)
lines(sx, predict(linreg, data.frame(x = sx)))

Indometh nonlinear fit

Indo.1 <- Indometh[Indometh$Subject == 1, ]
with(Indo.1, plot(time,conc))
indometh.nls <- nls(conc ~ c0 * exp(- k * time),  # model
    data = Indo.1,                # dataset
    start = list(c0 = 1, k = 0))  # initial guess for the parameters
lines(seq(0,8, by = 0.1),predict(indometh.nls, list(time = seq(0,8, by = 0.1))))

Nonlinear regression with an ODE model

rhs <- function(Time, conc, k) {
    return(list(c(- k * conc)))
}

odemodel <- function(Time,c0,k) {
  timespan <- sort(Time)
  addedzero <- FALSE
  if (timespan[1] > 0)
  {
    timespan <- c(0, timespan)
    addedzero <- TRUE
  }
  expmodel.sol <- ode(c(conc = c0), timespan, rhs, k)
  if (addedzero)
  {
    return( expmodel.sol[-1,2])
  }
  else
  {
    return(expmodel.sol[,2])
  }
}

with(Indo.1, plot(time,conc))
indometh.ode <- nls(conc ~ odemodel(time, c0, k),  # model
    data = Indo.1,                  # dataset
    start = list(c0 = 5, k = 1.2))  # initial guess for the parameters

Refining and selecting a model

biexpmodel <- function(Time,c01,c02,k1,k2) {
    return(c01*exp(- k1 * Time) + c02*exp(- k2 * Time))
}
with(Indo.1, plot(time,conc))
indometh.biexp <- nls(conc ~ biexpmodel(time,c01,c02,k1,k2),  # model
    data = Indo.1,                # dataset
    start = list(c01 = 2, c02 = 1, k1 = 2, k2 = 1))  # initial guess for the parameters
lines(seq(0,8, by = 0.1),predict(indometh.biexp, list(time = seq(0,8, by = 0.1))))
print(c(AIC_biexp = AIC(indometh.biexp), AIC_exp = AIC(indometh.nls)))
LS0tCnRpdGxlOiAiSW5mZXJyaW5nIG1vZGVsIHBhcmFtZXRlcnMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMgUG9seW5vbWlhbCBpbnRlcnBvbGF0aW9uCgpHZW5lcmF0aW9uIG9mIHRoZSByYXcgZGF0YQoKYGBge3J9Cm4gPC0gMTUKc2V0LnNlZWQoNzY1NDMyMSkKeCA8LSBydW5pZihuKSAjIGdlbmVyYXRlIG4gcmFuZG9tIG51bWJlcnMgaW4gWzAsMV0KeSA8LSBydW5pZihuLCBtaW4gPSAtMSwgMSkgIyBnZW5lcmF0ZSByZXNwb25zZSBpbiBbLTEsMV0KcGxvdCh4LHkpCmBgYAoKUG9seW5vbWlhbCBpbnRlcnBvbGF0aW9uIG9mIHRoZSByYXcgZGF0YQoKYGBge3J9CmludGVycCA8LSBsbSh5IH4gcG9seSh4LG4tMSkpICMgcG9seW5vbWlhbCBpbnRlcnBvbGF0aW9uCnBsb3QoeCx5KQpzeCA8LSBzb3J0KHgpCmxpbmVzKHN4LCBwcmVkaWN0KGludGVycCwgZGF0YS5mcmFtZSh4ID0gc3gpKSkKYGBgCgojIyBQcmVkaWN0aW9uIG9mIGEgbmV3IGRhdGEgcG9pbnQKCkdlbmVyYXRlIGEgbmV3IGRhdGEgcG9pbnQKCmBgYHtyfQpuZXd4IDwtIHJ1bmlmKDEpIApuZXd5IDwtIHJ1bmlmKDEpICMgYWRkIGEgbmV3IGRhdGEgcG9pbnQgCmBgYAoKUHJlZGljdGlvbgoKYGBge3J9CnBsb3QoYygwLDEpLGMoLTEwLDEpLCB0eXBlID0gJ24nLCB4bGFiID0gJ3gnLCB5bGFiID0gJ3knKQpwb2ludHMoeCx5KQpwb2ludHMobmV3eCxuZXd5LCBwY2ggPSAxOSkKbGluZXMoc29ydChjKHgsbmV3eCkpLCBwcmVkaWN0KGludGVycCwgZGF0YS5mcmFtZSh4ID0gc29ydChjKHgsbmV3eCkpKSkpCnBvaW50cyhzb3J0KGMoeCxuZXd4KSksIHByZWRpY3QoaW50ZXJwLCBkYXRhLmZyYW1lKHggPSBzb3J0KGMoeCxuZXd4KSkpKSwgcGNoID0gMykKYGBgClRoZSBkaWZmZXJlbmNlIGJldHdlZW4gdGhlIGRhdGEgcG9pbnQgYW5kIHRoZSBwcmVkaWN0aW9uIGlzIGh1Z2UuCgpgYGB7cn0KcHJpbnQocHJlZGljdChpbnRlcnAsIGRhdGEuZnJhbWUoeCA9IG5ld3gpKSAtIG5ld3kpCmBgYAoKR3JhcGggb2YgdGhlIGZpdHRlZCBwb2x5bm9taWFsCgpgYGB7cn0KcGxvdChjKDAsMSksYygtMTAwLDEwMCksIHR5cGUgPSAnbicsIHhsYWIgPSAneCcsIHlsYWIgPSAneScpCnBvaW50cyh4LHkpCmxpbmVzKHNlcSgwLDEsIGJ5ID0gMC4wMDUpLCBwcmVkaWN0KGludGVycCwgZGF0YS5mcmFtZSh4ID0gc2VxKDAsMSwgYnkgPSAwLjAwNSkpKSkKYGBgCiMgTGluZWFyIHJlZ3Jlc3Npb24gb24gdGhlIHNhbWUgZGF0YXNldAoKYGBge3J9CmxpbnJlZyA8LSBsbSh5IH4geCkgIyBsaW5lYXIgcmVncmVzc2lvbgpwbG90KHgseSkKcG9pbnRzKG5ld3gsbmV3eSwgcGNoID0gMTkpCnN4IDwtIHNvcnQoeCkKbGluZXMoc3gsIHByZWRpY3QobGlucmVnLCBkYXRhLmZyYW1lKHggPSBzeCkpKQpgYGAKCiMgSW5kb21ldGggbm9ubGluZWFyIGZpdAoKYGBge3J9CkluZG8uMSA8LSBJbmRvbWV0aFtJbmRvbWV0aCRTdWJqZWN0ID09IDEsIF0Kd2l0aChJbmRvLjEsIHBsb3QodGltZSxjb25jKSkKaW5kb21ldGgubmxzIDwtIG5scyhjb25jIH4gYzAgKiBleHAoLSBrICogdGltZSksICAjIG1vZGVsCiAgICBkYXRhID0gSW5kby4xLCAgICAgICAgICAgICAgICAjIGRhdGFzZXQKICAgIHN0YXJ0ID0gbGlzdChjMCA9IDEsIGsgPSAwKSkgICMgaW5pdGlhbCBndWVzcyBmb3IgdGhlIHBhcmFtZXRlcnMKbGluZXMoc2VxKDAsOCwgYnkgPSAwLjEpLHByZWRpY3QoaW5kb21ldGgubmxzLCBsaXN0KHRpbWUgPSBzZXEoMCw4LCBieSA9IDAuMSkpKSkKYGBgCgojIE5vbmxpbmVhciByZWdyZXNzaW9uIHdpdGggYW4gT0RFIG1vZGVsCgpgYGB7cn0KcmhzIDwtIGZ1bmN0aW9uKFRpbWUsIGNvbmMsIGspIHsKICAgIHJldHVybihsaXN0KGMoLSBrICogY29uYykpKQp9CgpvZGVtb2RlbCA8LSBmdW5jdGlvbihUaW1lLGMwLGspIHsKICB0aW1lc3BhbiA8LSBzb3J0KFRpbWUpCiAgYWRkZWR6ZXJvIDwtIEZBTFNFCiAgaWYgKHRpbWVzcGFuWzFdID4gMCkKICB7CiAgICB0aW1lc3BhbiA8LSBjKDAsIHRpbWVzcGFuKQogICAgYWRkZWR6ZXJvIDwtIFRSVUUKICB9CiAgZXhwbW9kZWwuc29sIDwtIG9kZShjKGNvbmMgPSBjMCksIHRpbWVzcGFuLCByaHMsIGspCiAgaWYgKGFkZGVkemVybykKICB7CiAgICByZXR1cm4oIGV4cG1vZGVsLnNvbFstMSwyXSkKICB9CiAgZWxzZQogIHsKICAgIHJldHVybihleHBtb2RlbC5zb2xbLDJdKQogIH0KfQoKd2l0aChJbmRvLjEsIHBsb3QodGltZSxjb25jKSkKaW5kb21ldGgub2RlIDwtIG5scyhjb25jIH4gb2RlbW9kZWwodGltZSwgYzAsIGspLCAgIyBtb2RlbAogICAgZGF0YSA9IEluZG8uMSwgICAgICAgICAgICAgICAgICAjIGRhdGFzZXQKICAgIHN0YXJ0ID0gbGlzdChjMCA9IDUsIGsgPSAxLjIpKSAgIyBpbml0aWFsIGd1ZXNzIGZvciB0aGUgcGFyYW1ldGVycwpgYGAKCiMgUmVmaW5pbmcgYW5kIHNlbGVjdGluZyBhIG1vZGVsCgpgYGB7cn0KYmlleHBtb2RlbCA8LSBmdW5jdGlvbihUaW1lLGMwMSxjMDIsazEsazIpIHsKICAgIHJldHVybihjMDEqZXhwKC0gazEgKiBUaW1lKSArIGMwMipleHAoLSBrMiAqIFRpbWUpKQp9CndpdGgoSW5kby4xLCBwbG90KHRpbWUsY29uYykpCmluZG9tZXRoLmJpZXhwIDwtIG5scyhjb25jIH4gYmlleHBtb2RlbCh0aW1lLGMwMSxjMDIsazEsazIpLCAgIyBtb2RlbAogICAgZGF0YSA9IEluZG8uMSwgICAgICAgICAgICAgICAgIyBkYXRhc2V0CiAgICBzdGFydCA9IGxpc3QoYzAxID0gMiwgYzAyID0gMSwgazEgPSAyLCBrMiA9IDEpKSAgIyBpbml0aWFsIGd1ZXNzIGZvciB0aGUgcGFyYW1ldGVycwpsaW5lcyhzZXEoMCw4LCBieSA9IDAuMSkscHJlZGljdChpbmRvbWV0aC5iaWV4cCwgbGlzdCh0aW1lID0gc2VxKDAsOCwgYnkgPSAwLjEpKSkpCmBgYAoKYGBge3J9CnByaW50KGMoQUlDX2JpZXhwID0gQUlDKGluZG9tZXRoLmJpZXhwKSwgQUlDX2V4cCA9IEFJQyhpbmRvbWV0aC5ubHMpKSkKYGBg