Create a gist now

Instantly share code, notes, and snippets.

@dill /wiggly.R
Last active Oct 9, 2015

What would you like to do?
Goofy animation of integrating 2nd derivatives for a GAM penalty
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment