Goofy animation of integrating 2nd derivatives for a GAM penalty

wiggly.R

# David L Miller 2015, MIT license

library(numDeriv)
library(animation)

# function taken from the "Gu & Wahba 4 term additive example" from mgcv::gamSim
f2 <- function(x) 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10

xvals <- seq(0,1,len=100)

plot_wiggly <- function(f2, xvals){

  # pre-calculate
  f2v <- f2(xvals)
  f2vg <- grad(f2,xvals)
  f2vg2 <- unlist(lapply(xvals, hessian, func=f2))
  f2vg2min <- min(f2vg2) -2
  
  # now plot
  for(i in 1:length(xvals)){
    par(mfrow=c(1,3))
    plot(xvals, f2v, type="l", main="function", ylab="f")
    points(xvals[i], f2v[i], pch=19, col="red")
    
    plot(xvals, f2vg, type="l", main="derivative", ylab="df/dx")
    points(xvals[i], f2vg[i], pch=19, col="red")
    
    plot(xvals, f2vg2, type="l", main="2nd derivative", ylab="d2f/dx2")
    points(xvals[i], f2vg2[i], pch=19, col="red")
    polygon(x=c(0,xvals[1:i], xvals[i],f2vg2min),
            y=c(f2vg2min,f2vg2[1:i],f2vg2min,f2vg2min), col = "grey")
    
    ani.pause()
  }
}

saveGIF(plot_wiggly(f2, xvals), "wiggly.gif", interval = 0.2, ani.width = 800, ani.height = 400)