Hồi quy lượng tử (Làm tại nhà)


Vũ Sỹ Hoàng
1 năm trước
Hữu ích 0 Chia sẻ Viết bình luận 0
Đã xem 9173

Sau loạt bài viết của tôi  về các  thuật toán phân loại, đã đến lúc quay lại R, lần này là hồi quy lượng tử. Có, tôi vẫn muốn hiểu rõ hơn về thói quen tối ưu hóa trong R. Trước khi xem xét hồi quy lượng tử, chúng ta hãy tính toán trung vị hoặc lượng tử từ một mẫu.

Trung bình

Xét một mẫu  { y 1, ⋯, y n } . Để tính toán trung vị, giải:

có thể được giải quyết bằng các kỹ thuật lập trình tuyến tính. Chính xác hơn, vấn đề này tương đương với

một i , b i ≥0 và  y i - μ = một i - b i i = 1, ⋯, n .

Để minh họa, hãy xem xét một mẫu từ phân phối hợp lý,

n = 101 
set.seed(1)
y = rlnorm(n)
median(y)
[1] 1.077415


Đối với bài toán tối ưu hóa, hãy sử dụng dạng ma trận, với  3 n  ràng buộc và  2 n  tham số +1 :

library(lpSolve)
A1 = cbind(diag(2*n),0) 
A2 = cbind(diag(n), -diag(n), 1)
r = lp("min", c(rep(1,2*n),0),
rbind(A1, A2),c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y))
tail(r$solution,1) 
[1] 1.077415


Có vẻ như nó hoạt động tốt

Lượng tử

Tất nhiên, chúng ta có thể điều chỉnh mã trước đó cho các lượng tử:

tau = .3
quantile(x,tau)
      30% 
0.6741586


Chương trình tuyến tính bây giờ

với  một i , b i ≥0  và  y i - μ = một i - b ii = 1, ⋯, n . Mã R bây giờ là:

A1 = cbind(diag(2*n),0) 
A2 = cbind(diag(n), -diag(n), 1)
r = lp("min", c(rep(tau,n),rep(1-tau,n),0),
rbind(A1, A2),c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y))
tail(r$solution,1) 
[1] 0.6741586


Càng xa càng tốt…

Hồi quy lượng tử (Đơn giản)

Hãy xem xét các tập dữ liệu sau đây, với giá thuê căn hộ tại một thành phố lớn của Đức, như là một chức năng của bề mặt, năm xây dựng, v.v.

base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE)


Chương trình tuyến tính cho hồi quy lượng tử hiện nay

với  a i , b i ≥0  và

y i - [ β 0 τ + β 1 τ x i ] = một i - b i

i = 1, ⋯, n . Vì vậy, sử dụng ở đây

require(lpSolve) 
tau = .3
n=nrow(base)
X = cbind( 1, base$area)
y = base$rent_euro
A1 = cbind(diag(2*n), 0,0) 
A2 = cbind(diag(n), -diag(n), X) 
r = lp("min",
       c(rep(tau,n), rep(1-tau,n),0,0), rbind(A1, A2),
       c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) 
tail(r$solution,2)
[1] 148.946864   3.289674


Tất nhiên, chúng ta có thể sử dụng hàm R để phù hợp với mô hình đó:

library(quantreg)
rq(rent_euro~area, tau=tau, data=base)
Coefficients:
(Intercept)        area 
 148.946864    3.289674


Ở đây một lần nữa, nó có vẻ hoạt động khá tốt. Tất nhiên chúng ta có thể sử dụng một mức xác suất khác nhau và có được một âm mưu.

plot(base$area,base$rent_euro,xlab=expression(paste("surface (",m^2,")")),
     ylab="rent (euros/month)",col=rgb(0,0,1,.4),cex=.5)
sf=0:250
yr=r$solution[2*n+1]+r$solution[2*n+2]*sf
lines(sf,yr,lwd=2,col="blue")
tau = .9
r = lp("min",
       c(rep(tau,n), rep(1-tau,n),0,0), rbind(A1, A2),
       c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) 
tail(r$solution,2)
[1] 121.815505   7.865536
yr=r$solution[2*n+1]+r$solution[2*n+2]*sf
lines(sf,yr,lwd=2,col="blue")


Hồi quy lượng tử (bội số)

Bây giờ chúng ta đã hiểu làm thế nào để chạy chương trình tối ưu hóa với một hiệp, tại sao không thử với hai? Ví dụ, chúng ta hãy xem liệu chúng ta có thể giải thích tiền thuê căn hộ là hàm (tuyến tính) của bề mặt và tuổi của tòa nhà hay không.

require(lpSolve) 
tau = .3
n=nrow(base)
X = cbind( 1, base$area, base$yearc )
y = base$rent_euro
A1 = cbind(diag(2*n), 0,0,0) 
A2 = cbind(diag(n), -diag(n), X) 
r = lp("min",
       c(rep(tau,n), rep(1-tau,n),0,0,0), rbind(A1, A2),
       c(rep(">=", 2*n), rep("=", n)), c(rep(0,2*n), y)) 
tail(r$solution,3)
[1] 0.000000 3.257562 0.077501


Thật không may, lần này, nó không hoạt động tốt

library(quantreg)
rq(rent_euro~area+yearc, tau=tau, data=base)
Coefficients:
 (Intercept)         area        yearc 
-5542.503252     3.978135     2.887234


Kết quả khá khác nhau. Và, trên thực tế, một kỹ thuật khác có thể xác nhận cái sau ( IRLS  - Iterively Reweighted Least Squares)

eps = residuals(lm(rent_euro~area+yearc, data=base))
for(s in 1:500){
  reg = lm(rent_euro~area+yearc, data=base, weights=(tau*(eps>0)+(1-tau)*(eps<0))/abs(eps))
  eps = residuals(reg)
}
reg$coefficients
 (Intercept)         area        yearc 
-5484.443043     3.955134     2.857943


Tôi không thể tìm ra những gì đã xảy ra với chương trình tuyến tính. Không chỉ các hệ số rất khác nhau, mà các dự đoán cũng vậy:

yr = r$solution[2*n+1]+r$solution[2*n+2]*base$area+r$solution[2*n+3]*base$yearc
plot(predict(reg),yr)
abline(a=0,b=1,lty=2,col="red")



Bây giờ là lúc để điều tra trên mạng.

Hữu ích 0 Chia sẻ Viết bình luận 0
Đã xem 9173