2017-10-18 untitled from rstudio

untitled-2017-10-18.R

library(boot)
library(plyr)

set.seed(22042017) # for reproducibility
generate_data <- function(nsamples=100) {
  x <- rnorm(nsamples,mean=.75, sd=.5)
  y <-  4 + 5*x -3*x^2 + rnorm(length(x), mean=0, sd = 1)
  data.frame(x,y)
}


dat <- generate_data()
plot(dat, pch=20, col="grey")
curve(4 + 5*x -3*x^2, add=T, lwd=1.5)
abline(lm(y~x, data=dat), col="red", lwd=1.5)


# estimates from assumptions
summary(lm(y~x, data=dat))

# bootstrap estimates, no assumptions
bootfun <- function(data, index) {
  coef(lm(y~x, data=dat, subset = index))
}
boot(dat, bootfun, 10000)

# The "real" s.e. of beta estimates
sampling <- raply(10000, function() {
  dat <- generate_data()
  coef(lm(y~x, data=dat))
})
aaply(sampling, 2, sd)


# The residuals are clearly not normal
hist(resid(lm(y~x, data=dat)), prob=T, nclass=20)
sdev <- sd(resid(lm(y~x, data=dat)))
curve(dnorm(x, sd=sdev), add=T)