#
# Animation (simple version) of the conspiracy between nonlinearity and random X (A. Buja):
#
nonlinearity.randomX.animated <- function(slowness=200, new=F) {
    if(as.character(getRversion()) >= "2.14") on.exit(dev.flush(dev.hold(0)))
    mu <- function(x) x^2+1
    x.grid <- seq(0, 1, len=101)
    x.pr <- c(-.4,-.1)
    N <- 10
    x.new <- sort(runif(N))
    y.base <- 1.5
    y.0 <- c()
    y.1 <- c()
    if(new) dev.new(buffered=F, width=6, height=6)
    par(mgp=c(1,.5,0), mar=c(3,3,1,1))
    for(i in 1:500) {
        x.old <- x.new
        x.new <- sort(runif(N))
        for(a in seq(0,1,length=slowness)) {
            if(as.character(getRversion()) >= "2.14") dev.hold()
            x <- (1-a)*x.old + a*x.new
            y <- mu(x)
            m <- lm(y ~ x)
            f <- fitted(m)
            plot(x=c(-.5,1), y=c(0,2), type="n", xlab="x", ylab=bquote(y == bold(mu) ( x )), xaxt="n", yaxt="n")
            lines(x=x.grid, y=mu(x.grid), col="gray70", lwd=6)                 # gray curve
            abline(m, lwd=2)                                                   # fitted line
            abline(h=0, col="gray")                                            # hor axis
            lines(c(0,1), c(0,0), lwd=3, col="gray")                           # hor x range [0,1]
            abline(v=x.pr[1], col="gray")                                      # ver axis
            lines( x=c(rbind(x,x,NA)), y=c(rbind(pmax(f,y),0,NA)), col="gray") # ver lines to points
            points(x=x, y=y, pch=16, cex=1.2, col="blue")                      # blue points on curve
            points(x=x, y=rep(0,length(x)), pch=16, cex=1.2)                   # black points on x axis
            points(x=x, y=f, cex=1.5, col="red")                               # red points on fit
            if(length(y.0)>0) {
                points(x=rep(x.pr[1], length(y.0)), y=y.0, pch=3, col="red")
                lines(x=rep(c(x.pr[1],x.pr[2],NA), length(y.0)),
                      y=y.base+c(rbind(0,y.1-y.0,NA)), col="red")
            }
            mtext(text=paste("Dataset #",i,sep=""), side=1, line=1, at=par()$usr[1], adj=0.5)
            if(as.character(getRversion()) >= "2.14") dev.flush()
            ## Sys.sleep(.00)
        }
        if(as.character(getRversion()) >= "2.14") dev.hold()
        y.pr <- predict(m, newdata=data.frame(x=x.pr))
        y.0[i] <- y.pr[1]
        y.1[i] <- y.pr[2]
        lines(x=x.pr[c(1,2,2)], y=y.pr[c(1,1,2)], col="red", lwd=2)   # triangle lines
        points(x=x.pr, y=y.pr, col="red", cex=2, pch=16)              # 2 triangle points
        lines(x=x.pr, y=y.base+c(0,diff(y.pr)), col="black", lwd=2)        # triangle line from base
        points(x=x.pr, y=y.base+c(0,diff(y.pr)), col="red", cex=2, pch=16) # 2 triangle points from base
        if(as.character(getRversion()) >= "2.14") dev.flush()
        Sys.sleep(1.5)
    }
}
nonlinearity.randomX.animated()