###################
   x <- 1 : 8
   y <- x^2
   plot(x, y, xlim=c(0,9), ylim=c(0,70), type="l", lwd=2, col="blue", main="main_title", xlab="x_label", ylab="y_label")
###################

###################
   curve(exp(-x^(2)/2), from=-3, to=3, n=10, xlim=c(-3,3), ylim=c(0,1.2))
   curve(exp(-x^(2)/2), from=-3, to=3, n=100, xlim=c(-3,3), ylim=c(0,1.2))
###################

###################
   curve(exp(-x^(2)/2), from=-3, to=3, n=100, xlim=c(-3,3), ylim=c(0,1.2))
   par(new=T)
   curve(exp(-x^(2)/1), from=-3, to=3, n=100, xlim=c(-3,3), ylim=c(0,1.2), ann=F, col="blue", xlab="", ylab="")
   par(new=T)
   curve(exp(-x^(2)/0.5), from=-3, to=3, n=100, xlim=c(-3,3), ylim=c(0,1.2), ann=F, col="red", xlab="", ylab="")
   legend("topright", legend=c("var1", "var1/2", "var1/4"), lty=1, col=c("black", "blue", "red"))
###################

###################
   curve(exp(-x^(2)/2), from=-3, to=3, n=100, xlim=c(-3,3), ylim=c(0,1.2))
   abline(0,0)
   abline(v = 0)
###################

###################
   x <- seq(1,10, length=20)
   y <- x^2
   w <- x^(2.3)
   z <- x^(2.5)
   X <- cbind(y,w,z) ## matrix X
   matplot(x, X, type="o", pch=1:3, lty=1:3, col=1:3, xlab="", ylab="")
   legend("topleft", legend=c("x^2","x^2.3","x^2.5"), lty=1:3, pch=1:3, col=1:3)
###################

###################      
   x <- seq(1,10, length=20)
   y <- x^2
   v <- x^(2.1)
   w <- x^(2.3)
   z <- x^(2.5)
   par(mfrow=c(2,2))  ## let's try mfcol=c(2,2)
   plot(x,y,main=1); plot(x,v,main=2); plot(x,w,main=3); plot(x,z,main=4)  
###################   

   
   
###################      
   f <- function(x){exp(-abs(x))/2}
   curve(f, xlim=c(-3,3), ylim=c(0,0.5))
###################      

###################      
   rlap <- function(n){
   	R <- NULL; i <- 1; y <- 1/2
	while(i <= n){
		u1 <- runif(1,-10,10)  ## support is limitted
		u2 <- runif(1,0,y)
		if(u2 <= f(u1)){
			R <- c(R, u1); i <- i + 1
		}
	}
   return(R)
   }
###################         
   
   
###################   
   ns <- c(5, 10, 50, 100, 500, 1000, 5000, 10000)
   M <- matrix(0, 1000, length(ns))
   for(i in 1:length(ns)){
   	M[,i] <- means(ns[i], 1000) 
   }
   par(mfcol=c(1,4))
   boxplot(M, names=ns)
   hist((M[,1]-mean(M[,1]))/sd(M[,1]), main="", xlab="", ylab="")
   hist((M[,4]-mean(M[,4]))/sd(M[,4]), main="", xlab="", ylab="")
   hist((M[,8]-mean(M[,8]))/sd(M[,8]), main="", xlab="", ylab="")
###################   



###################   
   n <- dim(X)[1]; p <- dim(X)[2]
   set.seed(1); ran <- sample(1:n, 3000)
   trainX <- X[ran,]  
   testX <- X[-ran,]  ## trimming ran
   
   reg1 <- lm(formula = quality ~ volatile.acidity, data = trainX)
   summary(reg1)
   
   regful <- lm(formula = quality ~ ., data = trainX)
   summary(regful)
   
   y <- testX[,12]
   extendX <- cbind(rep(1,length(y)), testX[,-12])
   colnames(extendX)[1] <- "intercept"
   yhat1 <- as.matrix(extendX[, c(1,3)]) %*% as.vector(reg1$coef)
   yhatall <- as.matrix(extendX) %*% as.vector(regful$coef)
   mse1 <- sum((y-yhat1)^2)/length(y) ## mean squared error
   mseall <- sum((y-yhatall)^2)/length(y)
   
   table(y, round(yhat1))
   table(y, round(yhatall))
   
   ynew <- rep(0, length(n))
   yhat1new <- rep(0, length(n))
   yhatallnew <- rep(0,length(n)) 
   ynew[y >= 6] <- 1; ynew[y <= 5] <- 0
   yhat1new[round(yhat1) >= 6] <- 1; yhat1new[round(yhat1) <= 5] <- 0
   yhatallnew[round(yhatall) >= 6] <- 1; yhatallnew[round(yhatall) <= 5] <- 0
   table(ynew, yhat1new)
   table(ynew, yhatallnew)
   
   step(regful)               
###################   



###################   
cda <- function(y,x,init,lambda){
      n <- dim(x)[1]; p <- dim(x)[2]
      cov <- t(x)%*%x
      covy <- t(x)%*%y
      beta <- init
      judge <- 100
      z <- rep(0,p); old<-rep(0,p)
      repeat {
         if(judge >= 0.001){
	    for(j in 1:p){
	       old[j] <- beta[j]
	       z[j] <- covy[j] - sum(cov[j,]*beta) + cov[j,j]*beta[j]
	       beta[j] <- sign(z[j])*max(0,abs(z[j]) - lambda)/cov[j,j]
	     } 
	     judge <- sum(abs(beta-old))
	  } else break
      }
    return(beta)
}

x <- matrix(rnorm(100*20),100,20)
beta <- c(c(1,1,1),rep(0,17))
y <- x%*%beta + rnorm(100)
cda(y,x,rep(0,20),10)
###################   




###################   
pgda <- function(y,x,init,lambda,eta){
   n <- dim(x)[1]; p <- dim(x)[2]
   cov <- t(x)%*%x
   covy <- t(x)%*%y
   beta <- init
   judge <- 100
   repeat {
	if(judge >= 0.001){
	   old <- beta
	   bar_beta <- beta - (cov%*%beta - covy)/eta
	   beta <- sign(bar_beta)*pmax(rep(0,p), abs(bar_beta) - lambda/eta)
	   judge <- sum(abs(beta-old))
	   } else break
   }
   return(beta)	
}

x <- matrix(rnorm(100*20),100,20)
beta <- c(c(1,1,1),rep(0,17))
y <- x%*%beta + rnorm(100)
pgda(y,x,rep(0,20),10,100)
###################