### 4. CV i Bootstrap

library(ISLR2)
?Auto

n=dim(Auto)[1] ## 392
head(Auto)


### koristeci CV metodu procijenit cemo (ocekivanu) testnu gresku za procjenitelj dobiven linearnom regresijom
### model je E[mpg | x_1,...,x_p] = x%*%beta za x_k=horsepower^k za razne p>=1

plot(mpg~horsepower, data=Auto)

library(boot) ## sadrzi funkciju cv.glm
glm.fit <- glm(mpg ~ horsepower, data = Auto)  ## ako ne zadamo "family", glm radi standardnu linearnu regresiju

?cv.glm
cv.err <- cv.glm(Auto, glm.fit, K=10) ## po defaultu K=n, a funkcija gubitka ("cost") kvadratna greska
cv.err$delta  ## druga komponenta dobivena uz korekciju pristranosti

set.seed(17)
cv.error.10 <- rep(0, 10)
for (i in 1:10) {
  glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)
  cv.error.10[i] <- cv.glm(Auto, glm.fit, K = 10)$delta[1]
}
plot(1:10, cv.error.10, type="b")  
## nagli pad CV greske kada prelazimo na kvadratni model, ali onda ne izgleda da je potrebno polinome veceg
### bilo bi zanimljivo izracunati standardnu gresku i primijeniti one SE rule 
## (u ovom slucaju broj kovarijata predstavlja parametar kompleksnosnosti)


############## bootstrap ##################

## racunamo standardne greske za koeficijente u gornjem modelu

boot.fn <- function(data, index)
  {
  coef(lm(mpg ~ horsepower, data = data, subset = index))
  }
boot.fn(Auto, 1:392) ## koristimo sve podatke

## realizacija jednog bootstrap procjenitelja:
b=sample(392, size = 392, replace = TRUE)  ## dobijemo niz uniformnih indeksa na 1:392 
## bootstrap uzorak je sada Auto[b,] -> neki elementi originalnog uzorka ponovit ce se vise puta,
## a odgovarajuci bootstrap procjenitelj je
boot.fn(Auto, b) ## ponovi vise puta

## bootstrap procjena standardne greske
B=1000  ## moze i vise 
boot_procj=matrix(NA, nrow=B, ncol=2)
set.seed(19)
for (i in 1:B) {
  b=sample(392, size = 392, replace = TRUE)
  boot_procj[i,]=boot.fn(Auto, b) ## tu spremamo i-ti bootstrap procjenitelj (beta0 u prvi stupac, beta1 u drugi)
}

SE_boot=apply(boot_procj, 2, sd)
SE_boot
## [1] 0.849772753 0.007371543

## ili preko funkcije boot:
boot(Auto, boot.fn, R = B)
## dobivamo malo drugacije rezultate zbog drugacijih bootstrap uzoraka (mozemo povecati B)

# koristeci standardne formule
lm.fit=lm(mpg ~ horsepower, data = Auto)
summary(lm.fit)
## procjene standardnih gresaka su 0.717499 i 0.006446 -- ima razlike

attach(Auto)
plot(horsepower, mpg)  
abline(lm.fit)  ## veza je daleko od linearne -> bootstrap vjerojatno daje bolje procjene!

## kvadratni clan
boot.fn <- function(data, index){
  coef(
    lm(mpg ~ horsepower + I(horsepower^2), 
       data = data, subset = index)
  )
}
set.seed(1)
boot(Auto, boot.fn, 1000)

lm.fit2=lm(mpg ~ horsepower + I(horsepower^2), data = Auto)
summary(lm.fit2)$coef ## bolje poklapanje

t=seq(50,220, by=0.1)
points(t, predict(lm.fit2, data.frame(horsepower=t), type="response"), type="l")