Instantly share code, notes, and snippets.

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

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)
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.