Thực hành code_ Tuần 6

Tên: Mai Huy

MSSV: 43.01.104.065

Số thứ tự: 08

1) The Validation Set Approad

# Load thư viện ISLR
library(ISLR)
# set.seed dùng để tái tạo những vector random giống nhau theo tương ứng với giá trị được đưa vào hàm seed
set.seed (1)

# set.seed dùng để tái tạo những vector random giống nhau theo tương ứng với giá trị được đưa vào hàm seed
set.seed (1)
# Tạo 1 tập train có 196 quan sát được rút trích ngẫu nhiên từ 392 quan sát ban đầu
train = sample(392 ,196)

# Fit mô hình hồi quy đơn giản để dự đoán biến đầu ra mpg trong tập train
lm.fit=lm(mpg~horsepower,data=Auto,subset=train)

# attach dùng để khiến cho những biến feature trong dữ liệu có sẵn trong Rstudio theo tên
attach(Auto)

The following objects are masked from Auto (pos = 3):

    acceleration, cylinders, displacement, horsepower, mpg, name, origin, weight, year

The following objects are masked from Auto (pos = 4):

    acceleration, cylinders, displacement, horsepower, mpg, name, origin, weight, year

The following objects are masked from Auto (pos = 5):

    acceleration, cylinders, displacement, horsepower, mpg, name, origin, weight, year

# Sử dụng hàm predict để ước lượng giá trị output cho 392 quan sát và sử dụng hàm mean() để tính MSE của 196 quan sát trong tập validation
mean((mpg -predict (lm.fit ,Auto))[-train ]^2)

[1] 23.26601

# Fit mô hình hồi quy bậc 2 để dự đoán biến đầu ra mpg trong tập train
lm.fit2=lm(mpg∼poly(horsepower ,2),data=Auto , subset=train)

# Tính giá trị MSE cho tập validation gồm 196 quan sát
mean((mpg -predict (lm.fit2 ,Auto ))[- train]^2)

[1] 18.71646

# Fit mô hình hồi quy bậc 3 để dự đoán biến đầu ra mpg trong tập train
lm.fit3=lm(mpg∼poly(horsepower ,3),data=Auto , subset=train)

# Tính giá trị MSE cho tập validation gồm 196 quan sát
mean((mpg -predict (lm.fit3 ,Auto ))[- train]^2)

[1] 18.79401

Validation set error cho mô hình hàm bậc 2 và 3 lần lượt là 18.71 và 18.79. Tuy nhiên khi chúng ta có tập training ngẫu nhiên khác, những con số này cũng sẽ thay đổi

# set.seed dùng để tái tạo những vector random giống nhau theo tương ứng với giá trị được đưa vào hàm seed
set.seed(2)
# Tạo 1 tập train có 196 quan sát được rút trích ngẫu nhiên từ 392 quan sát ban đầu
train = sample(392 ,196)

# Fit mô hình hồi quy bậc 1 để dự đoán biến đầu ra mpg trong tập train
lm.fit=lm(mpg∼horsepower ,subset=train)

# Tính giá trị MSE cho tập validation gồm 196 quan sát
mean((mpg -predict (lm.fit ,Auto))[-train ]^2)

[1] 25.72651

# Fit mô hình hồi quy bậc 2 để dự đoán biến đầu ra mpg trong tập train
lm.fit2=lm(mpg∼poly(horsepower ,2),data=Auto , subset=train)

# Tính giá trị MSE cho tập validation gồm 196 quan sát
mean((mpg -predict (lm.fit2 ,Auto ))[- train]^2)

[1] 20.43036

# Fit mô hình hồi quy bậc 3 để dự đoán biến đầu ra mpg trong tập train
 lm.fit3=lm(mpg∼poly(horsepower ,3),data=Auto , subset=train)

# Tính giá trị MSE cho tập validation gồm 196 quan sát
mean((mpg -predict (lm.fit3 ,Auto ))[- train]^2)

[1] 20.38533

Validation set error cho mô hình hồi quy bậc 1,2,3 lần lượt là 25,72, 20.43, 20.38. Chúng ta thấy mô hình bậc 2 cải thiện sai số rất nhiều so với bậc 1 còn mô hình bậc 3 thì sai số cải thiện không quá đáng kể.

2) Leave-One-Out Cross-Validation

# Fit mô hình hồi quy tuyến tính để dự đoán biến đầu ra mpg khi biến đầu vào là horsepower trong tập train
glm.fit = glm(mpg~horsepower ,data= Auto)

# Tính các hệ số B0 và B1 của mô hình
coef(glm.fit)

(Intercept)  horsepower 
 39.9358610  -0.1578447 

Ngoài cách viết trên chúng ta có thể bằng hàm lm() quen thuộc

# Fit mô hình hồi quy tuyến tính để dự đoán biến đầu ra mpg khi biến đầu vào là horsepower trong tập train
 lm.fit=lm(mpg∼horsepower ,data=Auto)

# Tính các hệ số B0 và B1 của mô hình
coef(lm.fit)

(Intercept)  horsepower 
 39.9358610  -0.1578447 

# Load thư viện boot
library(boot)
# Fit mô hình hồi quy tuyến tính để dự đoán biến đầu ra mpg khi biến đầu vào là horsepower trong tập train
glm.fit=glm(mpg∼horsepower ,data=Auto)

# Hàm cv.glm() dùng để tính sai số K-fold cross-validation cho mô hình tuyến tính
cv.err=cv.glm(Auto ,glm.fit)

# Hiển thị 2 sai số, sai số đầu là sai số cross-validation, sai số thứ 2 là sai số hiệu chỉnh khi đền bù cho thiên vị (bias)
cv.err$delta

[1] 24.23151 24.23114

# Tạo ra 1 vector có 5 giá trị = 0
cv.error=rep(0,5)

# Hiển thị các sai số cross-validation của 5 mô hình từ bậc 1 -> 5
cv.error

[1] 24.23151 19.24821 19.33498 19.42443 19.03321

Không có sự cải thiện thực sự đáng kể khi tăng bậc của mô hình.

3) k-Fold Cross-Validation

Chúng ta sẽ fit từng mô hình từ bậc 1 cho tới bậc 10 và lưu trữ giá trị sai số cross-validaiton trong vector cv.error.10

# set.seed dùng để tái tạo những vector random giống nhau theo tương ứng với giá trị được đưa vào hàm seed
set.seed(17)
# Tạo ra 1 vector có 10 giá trị = 0
cv.error.10=rep(0 ,10)
# Với phương pháp k-FOld CV này, chúng ta thử đặt k=10 tức tập dữ liệu sẽ chia làm 10, 9 phục vụ cho tập train, và 1 cho tập validation, test error sẽ được tính bằng cách lấy trung bình của 10 MSE.
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]
}

# Hiển thị các sai số cross-validation của 10 mô hình từ bậc 1 -> 10
cv.error.10

 [1] 24.27207 19.26909 19.34805 19.29496 19.03198 18.89781 19.12061 19.14666 18.87013 20.95520

Khi dùng k-fold Cross-Validation thì tốc độ tính toán nhanh hơn nhiều so với LOOCV, và ở đây chúng ta Không có sự cải thiện thực sự đáng kể khi tăng bậc của mô hình khi sai số dao động tăng giảm liên tục.

4) The Bootstrap

Tính toán độ chính xác thống kê

# Tạo ra 1 hàm để tính toán độ chính xác thống kê (α) khi được truyền vào data và các quan sát được sử dụng cho việc tính toán α
alpha.fn=function(data,index){
X=data$X[index]
Y=data$Y[index]
return((var(Y)-cov(X,Y))/(var(X) + var(Y)-2*cov(X,Y)))
}

# Tính toán α sử dụng 100 quan sát đầu trong dữ liệu Portfolio
alpha.fn(Portfolio ,1:100)

[1] 0.5758321

# set.seed dùng để tái tạo những vector random giống nhau theo tương ứng với giá trị được đưa vào hàm seed
set.seed (1)
# Sử dụng hàm sample() để chọn 100 quan sát quan sát theo thứ tự ngẫu nhiên từ 100 quan sát ban đầu để tính toán α
alpha.fn(Portfolio,sample(100,100,replace=T))

[1] 0.7368375

Sử dụng hàm boot() để phân tích bootstrap

# Tạo ra 1000 boostrap để tính toán α
boot(Portfolio ,alpha.fn,R=1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Portfolio, statistic = alpha.fn, R = 1000)


Bootstrap Statistics :
     original       bias    std. error
t1* 0.5758321 -0.001695873  0.09366347

Hàm cho chúng ta thông tin về α = 0.5758 cho dữ liệu ban đầu, mức thiên vị (bias) và sai số chuẩn bootstrap SE(ˆα) = 0.09366.

Tính toán độ chính xác mô hình hồi quy

# Tạo ra hàm để tính 2 hệ số B0 và B1 của mô hình hồi quy tuyến tính dự đoán giá trị đầu ra là mpg và đầu vào là horsepower với số lượng quan sát tương ứng với index được truyền vào
boot.fn=function(data,index)
return(coef(lm(mpg~horsepower,data=data,subset=index)))

# Tính 2 hệ số B0 và B1 của mô hình hồi quy tuyến tính dự đoán giá trị đầu ra là mpg và đầu vào là horsepower với 392 quan sát đầu trong dữ liệu Auto
boot.fn(Auto ,1:392)

(Intercept)  horsepower 
 39.9358610  -0.1578447 

# set.seed dùng để tái tạo những vector random giống nhau theo tương ứng với giá trị được đưa vào hàm seed
set.seed (1)
# Sử dụng hàm sample() để chọn 392 quan sát quan sát theo thứ tự ngẫu nhiên từ 392 quan sát ban đầu để tính toán B0 và B1
boot.fn(Auto,sample(392,392,replace=T))

(Intercept)  horsepower 
 40.3404517  -0.1634868 

# Sử dụng hàm sample() để chọn 392 quan sát quan sát theo thứ tự ngẫu nhiên từ 392 quan sát ban đầu để tính toán B0 và B1
boot.fn(Auto,sample(392,392,replace=T))

(Intercept)  horsepower 
 40.1186906  -0.1577063 

# Sử dụng hàm boot() để tính toán sai số chuẩn của 1000 bootstrap cho 2 hệ số B0 và B1
boot(Auto,boot.fn,1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
      original        bias    std. error
t1* 39.9358610  0.0544513229 0.841289790
t2* -0.1578447 -0.0006170901 0.007343073

Chúng ta thấy rằng sai số chuẩn của B0 và B1 lần lượt là : SE(βˆ0) = 0.84129, SE(βˆ1) is 0.00734.

# Phân tich các hệ số của mô hình tuyến tính ở trên
summary(lm(mpg~horsepower,data=Auto))$coef

              Estimate  Std. Error   t value      Pr(>|t|)
(Intercept) 39.9358610 0.717498656  55.65984 1.220362e-187
horsepower  -0.1578447 0.006445501 -24.48914  7.031989e-81

Giá trị dự đoán mpg = 39.9358610 - 0.1578447 x horsepower, p-value của horsepower tương đối lớn nên không có bằng chứng thống kê cho thấy có mối quan hệ giữa biến horsepower và mpg.

# Tạo hàm tính sai số chuẩn boostrap của mô hình bậc 2 dự đoán giá trị đầu ra là mpg và đầu vào là horsepower với số lượng quan sát tương ứng với index được truyền vào
boot.fn=function (data ,index)
coefficients(lm(mpg∼horsepower +I(horsepower ^2),data=data ,
subset=index))

# set.seed dùng để tái tạo những vector random giống nhau theo tương ứng với giá trị được đưa vào hàm seed
set.seed(1)
# Sử dụng hàm boot() để tính toán sai số chuẩn của 1000 bootstrap cho 3 hệ số B0,B1, B3 
boot(Auto,boot.fn ,1000)

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
        original        bias     std. error
t1* 56.900099702  3.511640e-02 2.0300222526
t2* -0.466189630 -7.080834e-04 0.0324241984
t3*  0.001230536  2.840324e-06 0.0001172164

Chúng ta thấy rằng sai số chuẩn boostrap của B0, B1, B3 lần lượt là : SE(βˆ0) = 2.0300222526, SE(βˆ1) = 0.0324241984, SE(βˆ2) = 0.0001172164.