## ----include=FALSE------------------------------------------------------------
##### _Advanced Data Analysis from an Elementary Point of View_ #####
# R code for the chapter "Optimization"
# Please do not re-distribute or use without attribution
# http://www.stat.cmu.edu/~cshalizi/ADAfaEPoV/




## -----------------------------------------------------------------------------
my.newton = function(f,f.prime,f.prime2,beta0,tolerance=1e-3,max.iter=50) {
  beta = beta0
  old.f = f(beta)
  iterations = 0
  made.changes = TRUE
  while(made.changes & (iterations < max.iter)) {
   iterations <- iterations +1
   made.changes <- FALSE
   new.beta = beta - f.prime(beta)/f.prime2(beta)
   new.f = f(new.beta)
   relative.change = abs(new.f - old.f)/old.f -1
   made.changes = (relative.changes > tolerance)
   beta = new.beta
   old.f = new.f
  }
  if (made.changes) {
    warning("Newton's method terminated before convergence")
  }
  return(list(minimum=beta,value=f(beta),deriv=f.prime(beta),
              deriv2=f.prime2(beta),iterations=iterations,
              converged=!made.changes))
}


## ----regression-through-origin-scatterplot,echo=FALSE-------------------------
x <- runif(n=100,min=-1,max=1)
beta.true <- 4
y <- beta.true*x + rt(n=100,df=2)
plot(y~x)
abline(0,beta.true,col="grey")
abline(lm(y~x),lty=2)



## ----mse-curve,echo=FALSE,out.width="0.5\\textwidth"--------------------------
demo.mse <- function(b) { return(mean((y-b*x)^2)) }
curve(Vectorize(demo.mse)(x),from=0,to=10,xlab=expression(beta),ylab="MSE")
rug(x=beta.true,side=1,col="grey")



## ----lambda-vs-constraint,echo=FALSE,width="0.5\\textwidth"-------------------
lambda.from.c <- function(c) { mean(x*y)/sqrt(c) - mean(x^2) }
curve(lambda.from.c(x),from=0,to=20,xlab="c",ylab=expression(lambda))
abline(h=0, lty="dotted")



## ----penalized-beta-and-beta-squared, echo=FALSE------------------------------
par(mfrow=c(2,1))
beta.from.lambda <- function(l) { return(mean(x*y)/(l+mean(x^2))) }
curve(beta.from.lambda(x),from=0,to=6,
  xlab=expression(lambda),ylab=expression(tilde(beta)))
curve(beta.from.lambda(x)^2,from=0,to=6,
  xlab=expression(lambda),ylab=expression(tilde(beta)^2))
par(mfrow=c(1,1))