Resampling Methods
Pham Dinh Khanh
April 22, 2018
1 Resampling Methods
Resampling dữ liệu là một trong những phương pháp cần thiết trong thống kê hiện đại. Đó là việc chúng ta lặp đi lặp lại quá trình rút mẫu từ tập tập data tổng thể và hồi qui lại model để khảo sát tính hợp lý của model trong những điều kiện dữ liệu khác nhau. Bởi nếu chỉ hồi qui model trên một tập dữ liệu và đưa ra kết luận thì thường model sẽ không thể cover được hết tất cả các khả năng và không mang tính đại diện. Vì vậy việc lập đi lập lại model trên những tập dữ liệu được resampling sẽ đảm bảo tính chính xác và đại diện của model cao hơn và hạn chế được các trường hợp thông tin bị mất mát từ những quan sát không bao gồm trong training dataset nếu chỉ hồi qui 1 lần.
Có 2 phương pháp resampling chính được sử dụng trong machine learning bao gồm: cross validation và bootstrap. Cả 2 phương pháp đều quan trọng. Chẳng hạn như cross validation có thể được sử dụng để ước lượng sai số trong việc lựa chọn nhiều model với độ linh hoạt khác nhau. Còn bootstrap được sử dụng trong một vài tình huống, chẳng hạn như cung cấp một giải pháp để đánh giá mức độ chính xác của model cho một phương pháp machine learning nào đó.
2 Cross-validation
Như chúng ta đã biết, trong model hồi qui tuyến tính chúng ta thường sử dụng MSE (mean square error) hay còn gọi là bình quân sai số ngẫu nhiên để đo lường mức độ phụ hợp của model. Tuy nhiên MSE thường khác nhau đối với những tập training khác nhau nên nếu chúng ta chỉ lựa chọn hồi qui 1 lần trên một tập training và đưa ra kết luận về tính phù hợp của model sẽ dẫn tới mắc sai lầm.
Xét mối quan hệ giữa mpg và horsepower theo kiểu tuyến tính và kiểu phương trình bậc 2.
library(ggplot2)
library(ISLR)
ggplot(Auto,aes(horsepower,mpg)) +
geom_point(colour="blue", shape = 21, size = 2, fill = "white", stroke = 1) +
geom_smooth(method='lm',formula = y ~ poly(x,1), level = 0, size = 0.5, aes(colour = "Linear")) +
geom_smooth(method='lm',formula = y ~ poly(x,2), level = 0, size = 0.5, aes(colour = "Degree 2")) +
scale_colour_manual(name="Poly type", values=c("red", "green")) +
labs(
x = "Horsepower",
y = "Miles per gallon"
) +
theme_bw()
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'ymax'
## Note: no visible binding for global variable 'ymin'Ta có thể thấy mối quan hệ giữa horsepower và mpg theo phương trình bậc 2 là phù hợp hơn so với phương trình tuyến tính. . Chúng ta sẽ kiểm chứng điều này bằng cách chia mẫu thành 2 phần train và test bằng nhau, mỗi phần chiếm 196 quan sát. Hồi qui model trên tập train và từ đó tính toán MSE trên tập test.
library(dplyr)
#Viết phương trình tính MSE
MSE <- function(x,y){
sum((x-y)^2)/length(x)
}
df <- data.frame(MSE.Linear = numeric(), MSE.Quadratic = numeric())
set.seed(1)
for (i in 1:10000){
#Chia mẫu train và test ngẫu nhiên theo tỷ lệ 50:50
train <- Auto %>% sample_n(nrow(Auto)*0.5)
test <- setdiff(Auto,train)
#Hồi qui model linear theo tập train
lm <- lm(data=train,mpg ~ horsepower)
#Dự báo trên tập test
test_hat <- predict(lm, newdata = test)
#Tính MSE_linear
MSE_linear <- MSE(test$mpg, test_hat)
#Hồi qui model quadratic theo tập train
lm <- lm(data=train,mpg ~ poly(horsepower,2,raw = TRUE))
#Dự báo trên tập test
test_hat <- predict(lm, newdata = test)
#Tính MSE_quadratic
MSE_quadratic <- MSE(test$mpg, test_hat)
c <- c(MSE_linear, MSE_quadratic)
#Lưu các MSE vào df
df <- rbind(df,c)
}
names(df) <- c("MSE.Linear","MSE.Quadratic")
head(df)Đồ thị MSE của 2 model Linear và Quadratic:
#Thêm số thư tự cho df
df <- cbind(Order = 1:10000,df)
ggplot(df, aes(Order)) +
geom_line(aes(y=MSE.Linear,colour="MSE.Linear")) +
geom_line(aes(y=MSE.Quadratic,colour="MSE.Quadratic")) +
scale_colour_manual(name="MSE Type",values=c("brown","blue")) +
labs(
title="10000 MSE of model sampling according to Linear and Quadratic",
x = "Sample Order",
y = "MSE"
) +
geom_point(data = filter(df,MSE.Linear < MSE.Quadratic), aes(x=Order,y=MSE.Linear), shape=16, size=3, fill="black")
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'Ta thấy mặc dù có 9996/10000 các trường hợp đưa ra kết quả MSE.Linear > MSE.Quadratic tuy nhiên vẫn có 4 điểm hình tròn màu đen trên đồ thị có MSE.Linear < MSE.Quadratic. Đó là các điểm:
df %>% filter(MSE.Linear < MSE.Quadratic)Như vậy trong trường hợp chúng ta resampling mẫu rơi vào các trường hợp thiểu số này và đưa ra kết quả MSE.Linear < MSE.Quadratic. Chúng ta kết luận rằng mối quan hệ tuyến tính là phù hợp hơn mối quan hệ theo đa thưc bậc 2 mặc dù trên thực tế bậc 2 phù hợp hơn. Như vậy để khẳng định một mối quan hệ chắc chắn giữa các biến chúng ta phải thực hiện phương pháp Cross Validation trên nhiều mẫu khác nhau. Như vậy kết luận của chúng ta mới đáng tin cậy và đủ vững. Tuy nhiên nhược điểm của Cross Validation đó là:
Chúng ta thường phải chịu chi phí thời gian cho việc hồi qui nhiều lần các mẫu khác nhau và tính toán MSE trên tập test. Đặc biệt đối với trường hợp mẫu kích thước càng lớn thì tổ hợp các mẫu khác nhau càng lớn và thời gian dành cho mỗi chu kì tính toán cũng lớn hơn.
Việc resampling phải xác định trước tỷ lệ chia mẫu train:test là bao nhiêu sẽ phù hợp. Có nhiều quan điểm khác nhau cho vấn đề này. Có người sử dụng mẫu 50:50 hoặc quan điểm khác cho rằng 70:30 hoặc 80:20 sẽ hợp lý hơn. Tuy nhiên trong mọi trường hợp thì số quan sát của train luôn phải lơn hơn hoặc bằng test bởi vì các phương pháp đều có xu hướng cho thấy model sẽ dự báo kém hơn khi được train trên mẫu ít quan sát hơn.
3 Phương pháp LOOCV (Leave one out cross validation)
Giống như tên gọi của nó, LOOCV sẽ rút ra một quan sát bất kì trong mẫu tổng thể kích thước n và hồi qui model trên n-1 quan sát còn lại và sử dụng model đó để dự báo và tính MSE cho quan sát được rút. Quá trình này lập lại n lần từ quan sát thứ nhất đến quan sát thứ n. MSE ở lần rút quan sát thứ i được tính như sau: \[MSE_{i} = (y_{i}-\hat y_{i})^2\] Quá trình LOOCV sẽ ước lượng bình quân MSE của n lần rút mẫu: \[CV_{(n)}=\frac{1}{n} \sum_{i=1}^{n} MSE_{i}\] Phương pháp LOOCV có 2 ưu điểm so với những cách tiếp cận khác:
Số lượng quan sát sử dụng để fit model là n-1 và gần như là toàn bộ số quan sát nên kết quả của model fit phản ánh chính xác bản chất mối quan hệ dữ liệu hơn so với model hồi qui từ một nửa số quan sát.
Kết quả từ phương pháp Cross-Validation thông thường sẽ rất ngẫu nhiên do số lần thực hiện chia mẫu là ngẫu nhiên nên rất khó để đánh giá dựa trên kết quả đưa ra từ các phương pháp này. Trái lại LOOCV chỉ đưa ra một kết quả duy nhất và không bị thay đổi khi thực hiện resample mẫu là \(CV_{n}\) nên sử dụng chỉ số này sẽ ổn định hơn.
Sử dụng phương pháp LOOCV để đánh giá lại tính phù hợp giữa 2 model Linear và Quadratic:
df_loocv <- data.frame(MSE.Linear = numeric(), MSE.Quadratic = numeric())
set.seed(2)
for (i in 1:nrow(Auto)){
#Chia mẫu train và test ngẫu nhiên theo tỷ lệ 50:50
train <- Auto[-i,]
test <- setdiff(Auto,train)
#Hồi qui model linear theo tập train
lm <- lm(data=train,mpg ~ horsepower)
#Dự báo trên tập test
test_hat <- predict(lm, newdata = test)
#Tính MSE_linear
MSE_linear <- MSE(test$mpg, test_hat)
#Hồi qui model quadratic theo tập train
lm <- lm(data=train,mpg ~ poly(horsepower,2,raw = TRUE))
#Dự báo trên tập test
test_hat <- predict(lm, newdata = test)
#Tính MSE_quadratic
MSE_quadratic <- MSE(test$mpg, test_hat)
c <- c(MSE_linear, MSE_quadratic)
#Lưu các MSE vào df
df_loocv <- rbind(df_loocv,c)
}
names(df_loocv) <- c("MSE.Linear","MSE.Quadratic")
head(df_loocv)Giá trị của \(CV_{n}\) của model linear và model Quadratic:
cv <- df_loocv %>% summarise_at(vars(matches("MSE")),mean)
colnames(cv) <- c("CV.Linear","CV.Quadratic")
cvNhư vậy căn cứ vào CV.Linear > CV.Quadratic ta kết luận rằng model dạng Quadratic phù hợp hơn dạng Linear.
3.1 K-fold Cross validation
Một phương pháp thay thế cho LOOCV là k-fold CV. Phương pháp này chia mẫu thành k phần không overlap có kích thước gần như bằng nhau (trường hợp không chia hết thì sẽ gần bằng) gọi là các folds. Chọn ngẫu nhiên 1 fold là tập test, k-1 folds còn lại sẽ được sử dụng để fit model. Dựa trên model được fit, MSE sau đó sẽ được tính cho tập dữ liệu test. Ta lập lại quá trình này k lần trên các folds khác nhau và thu được các chỉ số MSE lần lượt là \(MSE_{1}, MSE_{2},..., MSE_{k}\). Giá trị k-fold CV được ước lượng như sau: \[CV_{k} = \frac{1}{k} \sum_{i=1}^{k} MSE_{i}\]
Ta có thể dễ dàng nhận ra LOOCV là trường hợp đặc biệt của K-fold CV với k = n. Trong thực tế chúng ta thường chọn k nhỏ hơn nhằm mục đích đẩy nhanh tốc độ tính toán hơn. Chẳng hạn trong ví dụ này ta sẽ sử dụng k=5 hoặc k=10. Như vậy chúng ta sẽ chỉ phải hồi qui 5-10 lần, ít hơn rất nhiều so với phương pháp LOOCV khi phải hồi qui model n lần. Với các mẫu kích thước lớn (chẳng hạn > 1 triệu quan sát) thì chi phí thời gian bỏ ra là đáng kể.
Đến đây chúng ta sẽ thắc mắc rằng tại sao phải chia làm k mẫu không overlap có kích thước gần bằng nhau và thực hiện k lần tính MSE mà không resample ngẫu nhiên k lần và tính toán MSE.
Thực hiện chia mẫu 10 lần một cách ngẫu nhiên theo tỷ lệ train:test là 50:50 và tính MSE trên mỗi tập train đó.
df_random <- data.frame()
set.seed(2)
for(i in 1:10){
train <- sample_n(Auto, nrow(Auto)*0.5)
test <- setdiff(Auto,train)
for(j in 1:10) {
#Hồi qui model theo đa thức bậc j
lm <- lm(data = train,mpg ~ poly(horsepower,j,raw=TRUE))
#Dự báo trên tập test
test_hat <- predict(lm, newdata = test)
#Tính MSE
MSE_value <- MSE(test$mpg, test_hat)
c <- c(i,j,MSE_value)
#Lưu vào df bậc của đa thức và MSE của tập test.
df_random <- rbind(df_random,c)
}
}
colnames(df_random) <- c("ResampleTime","Degree","MSE")Vẽ biểu đồ các giá trị MSE cho 10 lần resample mẫu theo các đa thức có bậc tăng dần.
ggplot(df_random,aes(x=Degree, y=MSE, colour=as.factor(ResampleTime))) +
geom_line() +
labs(
title = "MSE following 10 randomely resampling time corresponding with polygonal degree",
colour = "Resample Time"
) +
scale_x_continuous(breaks=1:10) +
ylim(0,60)
Ta thấy khi resample 10 lần ngẫu nhiên thì các đường MSE rất phân tán và không thể sử dụng ngẫu nhiên một giá trị MSE từ các lần resample này để đưa ra kết luận về mức độ phù hợp của model.
Resample 10 lần các mẫu theo phương pháp k-fold CV.
#Hàm số tính CV của MSE cho mỗi mẫu
CV_kfold <- function(k){
cv_kfold <- data.frame()
remain <- Auto
#Chia mẫu thành k phần
for(i in k:1){
fold_test <- sample_frac(remain,1/i)
fold_train <- setdiff(Auto,fold_test)
remain <- setdiff(remain, fold_test)
#Tính MSE cho tập được rút theo các đa thức bậc 1-10
for(j in 1:10){
lm <- lm(data=fold_train, mpg ~ poly(horsepower,j,raw=TRUE))
test_hat <- predict(lm, newdata = fold_test)
MSE_value <- MSE(fold_test$mpg, test_hat)
c <- c(i,j,MSE_value)
cv_kfold <- rbind(cv_kfold,c)
}
}
colnames(cv_kfold) <- c("ResampleTime","Degree","MSE")
cv_kfold <- cv_kfold %>% group_by(Degree) %>% summarise(CV=mean(MSE))
cv_kfold
}
#10-folds cv 10 lần ngẫu nhiên.
df_kfold <- data.frame()
for(i in 1:10){
df_kfold <- rbind(data.frame(ResampleTime=i,CV_kfold(10)),df_kfold)
}## Note: no visible binding for global variable 'Degree'df_kfold %>% head()Biểu đồ biểu diễn các giá trị MSE của phương pháp k-fold CV theo bậc của đa thức tăng dần:
ggplot(df_kfold, aes(x=Degree,y=CV,colour=as.factor(ResampleTime)))+
geom_line() +
labs(
title = "CV following 10-fold CV corresponding with polygonal degree",
colour = "Resample Time"
) +
scale_x_continuous(breaks=1:10) +
ylim(0,60)
Ta có thể thấy giá trị MSE của 10 lần resample khác nhau không quá khác biệt kể cả với những dạng model có bậc đa thức khác nhau. Điều đó cho thấy phương pháp k-fold CV rất ổn định trong việc đánh giá các sai số MSE. Đó là lý do chúng ta không sample mẫu ngẫu nhiên k lần mà sẽ chia thành k tập con có kích thước gần bằng nhau và sử dụng CV để đánh giá mức độ phù hợp của model.
3.2 Sự đánh đổi giữa Bias-Variance của phương pháp k-fold Cross-validation
Phương pháp Cross-validation thông thường thường dẫn tới overestimate giá trị MSE của tập test vì model fit chỉ hồi qui dựa trên 50% tập dữ liệu. Cũng dựa trên quan điểm này thì LOOCV sẽ trả về kết quả của MSE là xấp xỉ không chệch vì model fit được tính trên n-1 quan sát và gần như cover toàn bộ tập dữ liệu. Phương pháp k-fold sẽ đưa ra kết quả MSE không chệch ở mức trung bình vì mẫu được sử dụng cho model fit chỉ là (k-1)n/k và không cover toàn bộ mẫu. Khi k càng lớn thì MSE sẽ càng không chệch và ngược lại. Điều đó cho thấy nếu chúng ta muốn xét trên khía cạnh không chệch thì LOOCV là một lựa chọn phù hợp hơn so với k-fold CV. Tuy nhiên độ chệch không phải là yếu tố duy nhất mà chúng ta quan tâm trong quá trình ước lượng mà bên cạnh đó khía cạnh phương sai cũng được xem xét. Và thực tế cho thấy rằng LOOCV luôn có phương sai lớn hơn so với k-fold CV. Lý do là bởi phương pháp LOOCV sử dụng n-1 điểm dữ liệu để fit model nên dường như toàn bộ tập dữ liệu đã được sử dụng. Khi n đủ lớn ta có thể coi n phương trình fit model là tương đương nhau. Các sai số được tính ra do đó cũng sẽ tương quan. Còn với phương pháp k-fold CV thì k folds là tách biệt nhau và model fit được hồi qui trên những mẫu không trùng nhau quá lớn nên mức độ tương quan của các sai số được tính ra trên tập fold test sẽ ít tương quan hơn. Một mẫu có mức độ tương quan lớn sẽ dẫn đến phương sai lớn hơn nên phương pháp LOOCV sẽ có kết quả phương sai lớn hơn k-fold CV.
Đồ thị MSE của phương pháp loocv đối với phương trình linear.
#Sai số chuẩn của các MSE:
sd(df_loocv$MSE.Linear)## [1] 36.84434ggplot(data = df_loocv, aes(x=1:392,y=MSE.Linear)) +
geom_point(shape=21,size=2,colour='blue',fill='white',stroke=1) +
labs(
title="MSE according to LOOCV methology with linear model",
x="ResampleTime",
y="MSE"
) +
ylim(0,400)
Viết hàm tính MSE cho từng tập fold trong phương pháp k-fold CV.
MSE_kfold <- function(k){
cv_kfold <- data.frame()
remain <- Auto
#Chia mẫu thành k phần
for(i in k:1){
fold_test <- sample_frac(remain,1/i)
fold_train <- setdiff(Auto,fold_test)
remain <- setdiff(remain, fold_test)
#Tính MSE cho tập được rút theo đa thức bậc 1
lm <- lm(data=fold_train, mpg ~ horsepower)
test_hat <- predict(lm, newdata = fold_test)
MSE_value <- MSE(fold_test$mpg, test_hat)
c <- c(i,MSE_value)
cv_kfold <- rbind(cv_kfold,c)
}
colnames(cv_kfold) <- c("ResampleTime","MSE")
cv_kfold
}
df_MSEKfold <- MSE_kfold(40)Đồ thị MSE của phương pháp k-fold CV đối với phương trình linear.
#Sai số chuẩn của các MSE:
sd(df_MSEKfold$MSE)## [1] 12.18005ggplot(df_MSEKfold, aes(x=ResampleTime,y=MSE)) +
geom_point(shape=21,size=2,colour='blue',fill='white',stroke=1) +
labs(
title="MSE according to k-fold CV methology with linear model",
x="ResampleTime",
y="MSE"
) +
ylim(0,400)
Như vậy luôn có sự đánh đổi giữa sự không chệch-phương sai liên quan đến việc lựa chọn k trong phương pháp k-fold cross validation. Kết quả thực nghiệm cũng chỉ ra rằng khi chọn k=5 hoặc k=10 thì MSE cũng không quá chệch và quá biến động nên thông thường chúng ta có thể chọn các chỉ số này để cross validation.
3.3 Các vấn đề của Cross-Validation trong model phân loại.
Đối với các bài toán hồi qui có các đại lượng cần dự báo là biến định lượng thì chúng ta sử dụng MSE để đánh giá chất lượng của một phương pháp. Tuy nhiên với các bài toán phân loại chúng ta sẽ sử dụng một chỉ số khác là tỷ lệ phân loại nhầm Err và cách áp dụng thuật toán cross validation cũng gần như tương tự như trên. Chẳng hạn đối với thuật toán LOOCV thì giá trị trả về là tỷ lệ phân nhóm sai LOOCV có công thức như sau: \[CV_{n}=\frac{1}{n} \sum_{i=1}^{n} Err_{i}\]
Trong đó \(Err_{i} = I(y_{i} \neq \hat y_{i})\) tức là xác xuất để dự báo sai kết qủa. Trong trường hợp k-fold CV và cross validation thông thường cũng được thiết lập tương tự. Ở ví dụng bên dưới chúng ta sẽ sử dụng model logistic để dự báo loài hoa có phải là virginica hay loài khác dựa trên Sepal.Width và Sepal.Length của tập dữ liệu iris. Căn cứ vào bayes decision boundary chúng ta sẽ phân loại các quan sát thành 2 nhóm đối lập là Virginica và Other. Mỗi một phương trình logistic sẽ đưa ra một bayes decisioin boundary khác nhau có thể là tuyến tính, curved, flexiable. Rõ ràng rằng model logistic sẽ không có một bayes decision boundary đủ linh hoạt nếu chỉ có dạng tuyến tính bậc nhất và khi đó. Tuy nhiên chúng ta hoàn toàn có thể chuyển các bayes decision boundary này thành những đường non-linear chỉ bằng cách thay đổi bậc của đa thức trong phương trình dự báo. Chẳng hạn như dạng quadratic như sau:
\[log(\frac{p}{1-p})= \beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{1}^2+\beta_{3}X_{2}+\beta_{4}X_{2}^2\]
Khi đó đường bayes decesion boundary của model sẽ có dạng như sau:
#Tạo tập dữ liệu iris2
iris2 <- iris %>% transmute(Width = Sepal.Width, Length = Sepal.Length, Species = as.factor(ifelse(Species == "virginica",1,0)))
#Tạo hàm vẽ decisionBoundary
library(stringr)
decisionBoundary <- function(data, formula, modelType="logit", title = ""){
vars <- str_split(formula,"\\~") %>% unlist()
var_res <- str_trim(vars[1])
if(str_detect(vars[2],"poly")){
var_list <- vars[2] %>% str_split("\\+") %>% unlist() %>% str_trim() %>% str_split(",")
var_x1 <- var_list[1] %>% unlist()
var_x1 <- var_x1[1] %>% str_split("\\(") %>% unlist()
var_x1 <- var_x1[2]
var_x2 <- var_list[2] %>% unlist()
var_x2 <- var_x2[1] %>% str_split("\\(") %>% unlist()
var_x2 <- var_x2[2]
var_pred <- c(var_x1, var_x2)
} else {
var_pred <- str_split(vars[2],"\\+") %>% unlist() %>% str_trim()
}
#tao newData
mix_x1 <- seq(min(data[,c(var_pred[1])]),max(data[,c(var_pred[1])]),length = 100)
mix_x2 <- seq(min(data[,c(var_pred[2])]),max(data[,c(var_pred[2])]),length = 100)
#mix giua 2 tap
newData <- expand.grid(mix_x1,mix_x2)
colnames(newData) <- c(var_pred[1],var_pred[2])
#du bao tren tap newData
if (modelType == "logit") {
logit.fit <- glm(formula, data = data, family = binomial)
exp.grid <- predict(logit.fit,newData,type="link")
} else if (modelType == "knn") {
exp.grid <- knn(data[,c(var_pred[1],var_pred[2])], newData,as.factor(unlist(data[,c(var_res)])), k=15, prob=TRUE)
}
exp_prob <- matrix(exp.grid,length(mix_x1),length(mix_x2))
contour(mix_x1,mix_x2, exp_prob,
levels=0.5, main=
title,axes = FALSE)
points(data, col=ifelse(data[,c(var_res)]==1, "coral", "cornflowerblue"))
points(newData, pch=".", cex=1.2, col=ifelse(exp_prob>0.5, "coral", "cornflowerblue"))
box()
}
par(mfrow=c(2,2),mar=c(rep(1,4)))
decisionBoundary(iris2,formula ="Species ~ Width + Length", title = "Degree=1")## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'knn'
## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'decisionBoundary(iris2,formula ="Species ~ poly(Width,2,raw=TRUE) + poly(Length,2,raw=TRUE)", title="Degree=2")
decisionBoundary(iris2,formula ="Species ~ poly(Width,3,raw=TRUE) + poly(Length,3,raw=TRUE)",
title="Degree=3")
decisionBoundary(iris2,formula ="Species ~ poly(Width,4,raw=TRUE) + poly(Length,4,raw=TRUE)",
title="Degree=4")
Trên toàn bộ tập dữ liệu thì tỷ lệ dự báo sai số là:
Err <- function(data, newdata, formula, degree){
vars <- str_split(formula,"\\~") %>% unlist()
var_res <- str_trim(vars[1])
if(str_detect(vars[2],"poly")){
var_list <- vars[2] %>% str_split("\\+") %>% unlist() %>% str_trim() %>% str_split(",")
var_x1 <- var_list[1] %>% unlist()
var_x1 <- var_x1[1] %>% str_split("\\(") %>% unlist()
var_x1 <- var_x1[2]
var_x2 <- var_list[2] %>% unlist()
var_x2 <- var_x2[1] %>% str_split("\\(") %>% unlist()
var_x2 <- var_x2[2]
var_pred <- c(var_x1, var_x2)
} else {
var_pred <- str_split(vars[2],"\\+") %>% unlist() %>% str_trim()
}
logit.fit <- glm(formula, data = data, family = binomial)
logit.pre <- predict(logit.fit, newdata = newdata, type="response")
logit.pre <- ifelse(logit.pre > 0.5,1,0)
df <- data.frame(newdata[,c(var_res)],logit.pre)
err <- 1-nrow(filter(df,df[,1]==df[,2]))/nrow(df)
err
}
Err(iris2,iris2,formula = "Species ~ poly(Width,1,raw=TRUE)+poly(Length,1,raw=TRUE)")## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'## [1] 0.1933333Err(iris2,iris2,formula = "Species ~ poly(Width,2,raw=TRUE)+poly(Length,2,raw=TRUE)")## [1] 0.1866667Err(iris2,iris2,formula = "Species ~ poly(Width,3,raw=TRUE)+poly(Length,3,raw=TRUE)")## [1] 0.1733333Err(iris2,iris2,formula = "Species ~ poly(Width,4,raw=TRUE)+poly(Length,4,raw=TRUE)")## [1] 0.16Ta có thể thấy khi bậc cao dần thì tỷ lệ dự báo sai số càng thấp. Bậc 4 có tỷ lệ dự báo chính xác lớn nhất khi chỉ 16% dự báo sai, bậc 1 dự báo sai số lớn nhất là 19.3%. Tuy nhiên đây mới chỉ là kết quả hồi qui model trên 1 mẫu và liệu có đáng tin cậy nếu ta kết luận rằng bậc 4 luôn dự báo chính xác hơn. Do đó chúng ta cần thự hiện k-fold CV cho những model này trên nhiều tập dữ liệu khác nhau để đưa ra kết luận tương tự như đã làm với các model dự báo thông thường và so sánh bằng tỷ lệ \(CV_{k}\)
Resample 10 lần các mẫu theo phương pháp k-fold CV.
#Hàm tạo formula theo bậc
gen_formula <- function(i, formula){
str_replace_all(formula,"1",as.character(i))
}
#Hàm số tính CV của Err cho mỗi mẫu
CV_kfold <- function(k,data,formula){
cv_kfold <- data.frame()
remain <- data
#Chia mẫu thành k phần
for(i in k:1){
fold_test <- sample_frac(remain,1/i)
fold_train <- setdiff(data,fold_test)
remain <- setdiff(remain, fold_test)
#Tính Err cho tập được rút theo các đa thức bậc 1-10
for(j in 1:10){
Err <- Err(fold_train,fold_test,formula = gen_formula(j,formula))
c <- c(i,j,Err)
cv_kfold <- rbind(cv_kfold,c)
}
}
colnames(cv_kfold) <- c("ResampleTime","Degree","Err")
cv_kfold <- cv_kfold %>% group_by(Degree) %>% summarise(CV=mean(Err))
cv_kfold
}
#10-folds cv 10 lần ngẫu nhiên.
df_kfold <- data.frame()
for(i in 1:10){
df_kfold <- rbind(data.frame(ResampleTime=i,CV_kfold(10,iris2,"Species ~ poly(Width,1,raw=TRUE)+poly(Length,1,raw=TRUE)")),df_kfold)
}## Note: no visible binding for global variable 'Degree'df_kfold %>% head()Biểu đồ biểu diễn các giá trị CV của phương pháp k-fold CV theo bậc của đa thức tăng dần:
df_kfold %>% group_by(Degree) %>% summarise(CV=mean(CV)) %>%
ggplot(aes(x=Degree,y=CV)) +
geom_point(shape=21, size=2, fill = 'white', stroke=1, colour='blue') +
geom_line() +
labs(
title = "CV following 10-fold CV"
) +
scale_x_continuous(breaks=1:10) +
ylim(0,0.5)
Biểu đồ này cho thấy bậc 2 mới là model cho kết quả phân loại có tỷ lệ sai lầm nhỏ nhất.Như vậy nếu không cross-validation thì chúng ta có thể đã mắc sai lầm khi kết luận rằng bậc 4 sẽ cho kết quả dự báo chuẩn hơn.
4 bootstrap Method
Bootstrap là một phương pháp ứng dụng phổ biến và mạnh mẽ trong thống kê. Nó giúp ta đánh giá về các hệ số ước lượng hoặc các phương pháp thống kê học không chắc chắn. Chẳng hạn như, bootstrap có thể được sử dụng để đánh giá sai số chuẩn của hệ số ước lượng trong các model hồi qui tuyến tính. Mặc dù các phương pháp hồi qui trên thực tế đều đã tính toán sẵn sai số chuẩn của hệ số ước lượng. Tuy nhiên điểm mạnh của bootstrap là ở chỗ đối với các sai số khó ước lượng mà các phần mềm hoặc R package không tính sẵn thì bootstrap có thể ứng dụng để kiểm tra mức độ biến động của các đại lượng này. Chẳng hạn như một ví dụ nhỏ sau đây:
Giả sử rằng chúng ta muốn đầu tư một lượng tiền fixed vào một danh mục đầu tư có yield return lần lượt là X và Y, trong đó X, Y là các biến ngẫu nhiên. Giả định rằng tỷ lệ tiền đầu tư vào X và Y lần lượt là \(\alpha\) và \(1-\alpha\). Bài toán đặt ra cho chúng ta là lựa chọn được \(\alpha\) sao cho volatile hoặc mức biến động hoặc phương sai của danh mục là tối đa. Hay nói cách khác chúng ta cần tối thiểu hóa \(Var(\alpha X+(1-\alpha Y))\). Giá trị \(\alpha\) được đưa ra cho để tối thiểu hóa rủi ro danh mục như sau:
\[\hat \alpha = \frac{\hat \sigma_{Y} ^2- \hat \sigma_{XY}}{\hat \sigma_{X}^2+\hat \sigma_{Y}^2-2\hat\sigma_{X}\hat\sigma_{Y}}\]
Trong đó \(\hat \sigma^2\) là giá trị ước lượng của các phương sai (variance) và \(\hat\sigma_{XY}\) là hiệp phương sai (covariance) của X và Y.
Chúng ta sẽ thực hiện mô phỏng trên 100 cặp các giá trị return của các tài sản X và Y và thay vào công thức trên để ước lượng giá trị \(\alpha\). Từ đó tìm ra phân phối của các giá trị này. Phương sai của chúng không quá lớn so với giá trị kì vọng của \(\alpha\) thì ta có thể kết luận rằng tỷ lệ đầu tư ở mức \(E(\alpha)\) là hợp lý để tối thiểu hóa rủi ro.
Sử dụng package VNDS để lấy chuỗi lợi suất 2 mã chứng khoán VIC và VND
library(VNDS)
#Lấy dữ liệu 2 chuỗi VIC và VND
X <- VNDS::getSymbols("VIC",
"2017-04-26",
"2018-04-25",
src = "CP68")## #VIC from 2017-04-26 to 2018-04-25 already clonedY <- VNDS::getSymbols("VND",
"2017-04-26",
"2018-04-25",
src = "CP68")## Note: no visible binding for global variable 'DATE'
## Note: no visible binding for global variable 'DATE'
## #VND from 2017-04-26 to 2018-04-25 already cloned#Tạo chuỗi logreturn
library(tidyquant)
X <- X %>% tq_transmute(select = CLOSE,
mutate_fun = periodReturn,
period = "daily", type = "log")
colnames(X) <- c("date","X")
Y <- Y %>% tq_transmute(select = CLOSE,
mutate_fun = periodReturn,
period = "daily", type = "log")
colnames(Y) <- c("date","Y")
Return <- left_join(X,Y,by="date")
#Loại ngày đầu tiên
Return <- Return %>% filter(date != "2017-04-26")Viết hàm bootstrap dữ liệu lợi suất. parameter của hàm gồm k là số lần bootstrap và data source. Kết quả trả về của hàm là một dataframe có STT bootstrap và giá trị \(\alpha\) tính theo công thức trên.
#Tạo hàm simulate data
simulate <- function(data,size) {
#Tạo index ngẫu nhiên kích thước bằng size để lấy mẫu
index <- round(runif(size,0,nrow(data)))
#Khởi tạo df_bootstrap
df_bootstrap <- data.frame()
for (j in index){
df_bootstrap <- rbind(data[j,],df_bootstrap)
}
df_bootstrap
}
#Tạo hàm tính variance và correlation
var_cor <- function(data){
data %>% summarise(var_X = var(X),
var_Y = var(Y),
cov_XY = cov(X,Y))
}
#Tạo hàm tính alpha
alpha_cal <- function(data){
Var <- var_cor(data)
Var %>% summarise(alpha =
(var_Y - cov_XY)/(var_X + var_Y - 2*cov_XY))
}
#Tạp hàm bootstraping cho 100 phiên
bootstrap <- function(data,k){
df_result <- data.frame()
for (i in 1:k){
#Simulate mẫu
df_bootstrap <- simulate(Return,100)
#Tính alpha theo công thức
alpha <- alpha_cal(df_bootstrap)
df_result <- rbind(df_result ,c(STT = i, alpha = alpha$alpha))
colnames(df_result) <- c("STT","alpha")
}
df_result
}Biểu diễn đồ thị của X,Y khi simulation 4 lần với kích thước là 100.
simu_1 <- simulate(Return,100)
simu_2 <- simulate(Return,100)
simu_3 <- simulate(Return,100)
simu_4 <- simulate(Return,100)
#Tạo hàm vẽ đồ thị của các return
plot_return <- function(data){
ggplot(data,aes(x=X,y=Y))+
geom_point(shape=21, size=2, stroke=1, fill='white', colour='blue') +
theme_bw()
}
par(mfrow=c(2,2),mar=c(rep(1,4)))
plot_return(simu_1)
plot_return(simu_2)
plot_return(simu_3)
plot_return(simu_4)
Ta thấy với 4 lần simulation thì biến động của X và Y cũng không quá lớn, giá trị kì vọng \(\alpha\) thu được từ 4 quá trình simulation này lần lượt là:
alpha_cal(simu_1)alpha_cal(simu_2)## Note: no visible binding for global variable 'var_Y'
## Note: no visible binding for global variable 'cov_XY'
## Note: no visible binding for global variable 'var_X'
## Note: no visible binding for global variable 'var_Y'
## Note: no visible binding for global variable 'cov_XY'alpha_cal(simu_3)alpha_cal(simu_4)alpha_cal(Return)Phương pháp booostrap sẽ khác so với cross-validation ở chỗ thay vì lấy lặp đi lặp lại những mẫu độc lập nhau từ tổng thể thì chúng ta sẽ tạo ra những mẫu phân biết nhau từ 1 mẫu gốc duy nhất. Cách làm này không đòi hỏi chúng ta phải tạo ra dữ liệu mới hay mở rộng tập dữ liệu mà từ một bộ dữ liệu gốc sẵn có đó ta có thể resample ra nhiều tập dữ liệu khác nhau và đánh giá sai số của model dựa trên các tập dữ liệu mới sinh này.
Chẳng hạn chúng ta có \(A\) là tập dữ liệu chỉ chứa 3 quan sát. Chúng ta chọn ngẫu nhiên n quan sát từ tập dữ liệu này nhằm mục đích tạo ra những mẫu bootstrap data set. Các mẫu này cho phép lặp lại các quan sát (hay còn gọi là lựa chọn replacement). Chẳng hạn như chúng ta có mẫu \(A^{1}\) lặp lại 2 lần quan sát thứ 3, 1 lần quan sát thứ 1 và không chứa quan sát thứ 2.
A <- data.frame(Obs=1:3,X=c(4.6,2.5,2.1),Y=c(3.9,6.5,1.7))
AA1 <- A[c(1,1,3),]
A1Dựa trên mẫu \(A^{1}\) ta sẽ ước lượng được giá trị cho \(\alpha\) gọi là \(\hat \alpha_{1}\). quá trình này lập lại n lần và ta sẽ thu được n mẫu bootstrap data sets là \(A^{1},A^{2},...,A^{n}\) và ứng với n ước lượng của \(\alpha\) gồm \(\hat\alpha_{1},\hat\alpha_{2},...,\hat\alpha_{n}\). Chúng ta có thể tính phương sai của các \(\alpha\) này như sau:
\[SE(\hat\alpha)=\frac{1}{n-1} \sqrt{\sum_{i=1}^{n}(\hat\alpha_{i}-\frac{1}{n}\sum_{i=1}^{n}\hat\alpha_{i})^2}\] Giá trị trên được coi như là ước lượng phương sai của \(\hat\alpha\) từ mẫu gốc. Tuy nhiên chúng ta không thể sử dụng giá trị này như ước phương sai cho \(\hat \alpha\) toàn bộ tổng thể vì từ mẫu gốc chúng ta không thể có thêm những dữ liệu mới từ tổng thể mà chỉ có các quan sát thu được trong phạm vi mẫu gốc này.
Bên dưới ta sẽ so sánh 2 phương pháp bootstrap thực hiện từ một mẫu gốc và từ toàn bộ tổng thể trên tập 1000 mẫu kích thước mỗi mẫu là 100 quan sát.
bootstrap trên tổng thể
#bootstrap 1000 lần từ mẫu gốc
set.seed(1)
bootstrap_alpha <- bootstrap(Return,1000)## Note: no visible global function definition for 'c'
## Note: no visible global function definition for 'c'#giá trị mean, variance alpha
mean_alpha <- mean(bootstrap_alpha$alpha)bootstrap trên một mẫu gốc
Tạo mẫu gốc từ tổng thể. Để mang tính đại diện ta sẽ lựa chọn mẫu sub_return sao cho giá trị \(\alpha\) trả về theo công thức trên xấp xỉ với mean_alpha của tổng thể (khoảng cách là 0.01).
for(i in 1:1000){
sub_return<- Return %>% sample_frac(0.5)
if(abs(alpha_cal(sub_return)-mean_alpha)<0.01){
break
} else {
print(i)
}
}## [1] 1
## [1] 2
## [1] 3#bootstrap 1000 lần từ mẫu gốc đơn vừa tạo
set.seed(2)
bootstrap_alpha_sub <- bootstrap(sub_return,1000)
mean_alpha_sub <- mean(bootstrap_alpha_sub$alpha)library(reshape)
#Vẽ biểu đồ
par(mfrow=c(1,3),mar=c(rep(1,4)))
bootstrap_alpha %>%
ggplot(aes(alpha)) +
geom_histogram(bins = 10, colour='black', fill='blue', alpha = 0.3) +
geom_vline(xintercept = mean_alpha, size = 1, colour='red') +
labs(
x="alpha",
y=""
) +
theme_classic()
bootstrap_alpha_sub %>%
ggplot(aes(alpha)) +
geom_histogram(bins = 10, colour='black', fill='green', alpha = 0.3) +
geom_vline(xintercept = mean_alpha_sub, size = 1, colour='red') +
labs(
x="Tổng thể",
y="alpha"
) +
theme_classic()
## Note: no visible binding for global variable 'xend'
## Note: no visible binding for global variable 'yend'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'y'left_join(bootstrap_alpha,bootstrap_alpha_sub, by="STT") %>%
melt(id=c("STT")) %>%
ggplot(aes(variable, value, fill=variable)) +
geom_boxplot(show.legend = FALSE) +
labs(
x="Bootstrap",
y="alpha"
) +
scale_x_discrete(breaks=c("alpha.x","alpha.y"),
labels=c("Tổng thể","Bootstrap"))
## Note: no visible binding for global variable 'weight'
## Note: no visible binding for global variable 'xmin'
## Note: no visible binding for global variable 'xmax'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'size'Ta có thể nhận thấy phân phối của \(\alpha\) thu được từ mẫu gốc bootstrap và tổng thể không khác nhau nhiều khi trung bình hội tụ quan mức 0.69 và các khoảng ngũ phân vị có giá trị tương đối bằng nhau trên biểu đồ boxplot. Trong khi đó phương sai của \(\alpha\) được tính từ bootstrap và tổng thể lần lượt là 0.0051 so với 0.0055 cho thấy về mặt giá trị thì \(\alpha\) biến động không đáng kể.
var(bootstrap_alpha$alpha)## [1] 0.005721211var(bootstrap_alpha_sub$alpha)## [1] 0.005872686Như vậy bootstrap khá đáng tin cậy trong việc đánh giá \(\alpha\) và chúng ta có thể tin cậy rằng tỷ lệ tiền mặt đầu tư vào danh mục X là \(\alpha\) ở mức 0.69 là một tỷ lệ tối ưu hóa rủi ro.