### 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")