library("ggplot2", lib.loc="/usr/lib/R/library")
library("boot", lib.loc="/usr/lib/R/library")

makePlots <- function(stddev){
  set.seed(1)
  x <- rnorm(100,mean=0,sd=1)
  h<-hist(x)
  xfit<-seq(min(x),max(x),length=100)
  yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
  yfit<-yfit*diff(h$mids[1:2])*length(x)
  lines(xfit,yfit,col="blue")

  eps <- rnorm(100,mean=0,sd=stddev)
  h<-hist(eps)
  xfit<-seq(min(eps),max(eps),length=100)
  yfit<-dnorm(xfit,mean=mean(eps),sd=sd(eps))
  yfit<-yfit*diff(h$mids[1:2])*length(eps)
  lines(xfit,yfit,col="blue")

  y <- -1+0.5*x + eps
  y

  plot(x,y,main="X vs Y Scatter Plot")

  model <- lm(y ~ x)
  print(summary(model))
  abline(model,col="red", lwd=10)

  abline(c(-1,0.5),col="blue", lwd=5)

  legend("bottomright",c("Least Squares","Population", "Polynomial"),lwd=c(2.5,2.5,2.5),col=c("red","blue","green"),lty=c(1,1,1))

  quad <- lm(y ~ poly(x,2,raw=TRUE))
  quad
  fun <- function(x) quad$coefficients[3]*x^2+quad$coefficients[2]*x+quad$coefficients[1]
  curve(fun,col="green",add=TRUE, lwd=5)

  rsq <- function(formula, data, indices){
    data <- data.frame(data)
    d <- data[indices,]
    fit <- lm(formula, data=d)
    return(summary(fit)$r.square)
  }

  data <- data.frame(x,y)

  results <- boot(data=data,statistic=rsq, R=1000, formula=y~x)
  print(results)
  plot(results)

  return(model)
}
makePlots(0.25)
makePlots(0.1)
makePlots(0.5)