-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathch5.R
More file actions
62 lines (47 loc) · 1.14 KB
/
ch5.R
File metadata and controls
62 lines (47 loc) · 1.14 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
require(ISLR)
require(boot)
?cv.glm
plot(mpg~horsepower,data=Auto)
## LOOCV
glm.fit=glm(mpg~horsepower, data=Auto)
cv.glm(Auto,glm.fit)$delta #pretty slow (doesnt use formula (5.2) on page 180)
##Lets write a simple function to use formula (5.2)
loocv=function(fit){
h=lm.influence(fit)$h
mean((residuals(fit)/(1-h))^2)
}
## Now we try it out
loocv(glm.fit)
cv.error=rep(0,5)
degree=1:5
for(d in degree){
glm.fit=glm(mpg~poly(horsepower,d), data=Auto)
cv.error[d]=loocv(glm.fit)
}
plot(degree,cv.error,type="b")
## 10-fold CV
cv.error10=rep(0,5)
for(d in degree){
glm.fit=glm(mpg~poly(horsepower,d), data=Auto)
cv.error10[d]=cv.glm(Auto,glm.fit,K=10)$delta[1]
}
lines(degree,cv.error10,type="b",col="red")
## Bootstrap
## Minimum risk investment - Section 5.2
alpha=function(x,y){
vx=var(x)
vy=var(y)
cxy=cov(x,y)
(vy-cxy)/(vx+vy-2*cxy)
}
alpha(Portfolio$X,Portfolio$Y)
## What is the standard error of alpha?
alpha.fn=function(data, index){
with(data[index,],alpha(X,Y))
}
alpha.fn(Portfolio,1:100)
set.seed(1)
alpha.fn (Portfolio,sample(1:100,100,replace=TRUE))
boot.out=boot(Portfolio,alpha.fn,R=1000)
boot.out
plot(boot.out)