source('http://myweb.uiowa.edu/pbreheny/7600/s16/notes/fun.R')

# Slide 5
lasso <- function(theta, l) l*abs(theta)
dLasso <- function(theta, l=1) rep(l,length(theta))
ridge <- function(theta, l=1) l*theta^2
dRidge <- function(theta, l=1) l*abs(theta)

res <- 99
x <- seq(-5, 5, len=res)
xx <- seq(0, 5,len=res)
col <- c("#FF4E37FF", "#008DFFFF")

par(mfrow=c(1,2), mar=c(5, 5, 0.5, 0.5), mgp=c(2, 1, 0), oma=c(0,0,2,0))
matplot(x, cbind(lasso(x, 1), ridge(x, 0.3)), type="l", lwd=3, lty=1, col=col,
        xlab=expression(beta), ylab=expression(P(beta)), xaxt="n", yaxt="n", bty="l")
axis(1, at=0, label=0)
matplot(xx, cbind(dLasso(xx, 1), dRidge(xx, 0.3)), type="l", lwd=3, lty=1, col=col,
        xlab=expression(abs(beta)), ylab=expression(P*"'"(abs(beta))), xaxt="n", yaxt="n", bty="l")
axis(1, at=0, label=0)
axis(2, at=c(0,1), label=c(0, expression(lambda)), las=1)
toplegend(legend=c("Lasso", "Ridge"), lwd=3, col=col)

# Slide 21
par(mgp=c(2, 1, 0), mar=c(5, 5, 3, 0.5))
col <- c("#FF4E37FF", "#008DFFFF")
xx <- seq(-2.5, 2.5, len=99)
plot(xx, xx, lwd=3, type="l", col="gray", xlab="z", ylab=expression(hat(beta)(z)), bty="n", las=1)
lines(xx, xx*(abs(xx) > 1), lwd=3, col=col[1])
lines(xx, (xx - sign(xx))*(abs(xx) > 1), lwd=3, col=col[2])
toplegend(legend=c("Hard", "Soft"), lwd=3, col=col)

# Slide 23
lam <- 0.5
n <- 10
sig <- 1
se <- sig/sqrt(n)

par(mar=c(5,5,1.5,0.5), mfrow=c(1,2))
ymax <- 1.5
b  <- 0
plot(0, 0, type="n", xlim=c(-2,2), ylim=c(0,ymax), las=1, xlab=expression(hat(beta)), ylab="Density", bty="n")
bb <- seq(-2,0,len=49)
polygon(c(bb, rev(bb)), c(dnorm(bb-lam, b, se), rep(0, length(bb))), col="gray80", border=NA)
bb <- seq(0,2,len=49)
polygon(c(bb, rev(bb)), c(dnorm(bb+lam, b, se), rep(0, length(bb))), col="gray80", border=NA)
p <- pnorm((lam-b)/se) - pnorm((-lam-b)/se)
lines(0, p*ymax, type="h", lwd=1.5)
mtext(expression(beta[0]==0))

par(mar=c(5,0.5,1.5,5))
b  <- 1
ymax <- 1.6
plot(0, 0, type="n", xlim=c(-2,2), ylim=c(0,ymax), las=1, xlab=expression(hat(beta)), ylab="", bty="n", yaxt="n")
bb <- seq(-2,0,len=49)
polygon(c(bb, rev(bb)), c(dnorm(bb-lam, b, se), rep(0, length(bb))), col="gray80", border=NA)
bb <- seq(0,2,len=49)
polygon(c(bb, rev(bb)), c(dnorm(bb+lam, b, se), rep(0, length(bb))), col="gray80", border=NA)
p <- pnorm((lam-b)/se) - pnorm((-lam-b)/se)
lines(0, p*ymax, type="h", lwd=1.5)
axis(4, at=c(0,ymax/2,ymax), label=c(0, 0.5, 1), las=1)
mtext("Probability", 4, line=3)
mtext(expression(beta[0]==1))