2011-12-10 30 views
9

Tôi đang cố gắng hoàn thiện phương pháp so sánh hồi quy và PCA, lấy cảm hứng từ blog Cerebral Mastication cũng đã được thảo luận từ một góc khác trên SO. Trước khi tôi quên, rất cám ơn JD Long và Josh Ulrich vì phần lớn cốt lõi của điều này. Tôi sẽ sử dụng điều này trong một khóa học học kỳ tiếp theo. Xin lỗi điều này là dài!So sánh trực quan của hồi quy & PCA

CẬP NHẬT: Tôi đã tìm thấy một cách tiếp cận khác gần như hoạt động (vui lòng khắc phục nếu bạn có thể!). Tôi đã đăng nó ở phía dưới. Một cách tiếp cận thông minh hơn và ngắn hơn nhiều so với tôi đã có thể nghĩ ra!

Tôi về cơ bản đã làm theo các sơ đồ trước đó cho đến một điểm: Tạo dữ liệu ngẫu nhiên, tìm ra đường phù hợp nhất, vẽ số dư. Điều này được thể hiện trong đoạn mã thứ hai bên dưới. Nhưng tôi cũng đào xung quanh và viết một số chức năng để vẽ các đường bình thường đến một đường thẳng qua một điểm ngẫu nhiên (các điểm dữ liệu trong trường hợp này). Tôi nghĩ rằng những công việc này tốt, và chúng được thể hiện trong First Code Chunk cùng với bằng chứng họ làm việc.

Bây giờ, đoạn mã thứ hai hiển thị toàn bộ điều đang hoạt động bằng cách sử dụng luồng giống như @JDLong và tôi đang thêm hình ảnh của ô vẽ. Dữ liệu màu đen, đỏ là hồi quy với số dư màu hồng, màu xanh là PC thứ nhất và màu xanh lam nhạt nên là chuẩn mực, nhưng rõ ràng là không. Các hàm trong First Code Chunk vẽ các normals này có vẻ tốt, nhưng có điều gì đó không đúng với phần trình diễn: Tôi nghĩ rằng tôi phải hiểu nhầm một cái gì đó hoặc truyền các giá trị sai. Các tiêu chuẩn của tôi đi ngang, có vẻ như là một đầu mối hữu ích (nhưng cho đến nay, không phải với tôi). Có ai thấy có gì sai ở đây không?

Cảm ơn, điều này đã được vexing tôi trong một thời gian ... Plot showing problem

Mã Đầu tiên Chunk (Chức năng để Vẽ normals và Proof Họ làm việc):

##### The functions below are based very loosely on the citation at the end 

pointOnLineNearPoint <- function(Px, Py, slope, intercept) { 
    # Px, Py is the point to test, can be a vector. 
    # slope, intercept is the line to check distance. 

    Ax <- Px-10*diff(range(Px)) 
    Bx <- Px+10*diff(range(Px)) 
    Ay <- Ax * slope + intercept 
    By <- Bx * slope + intercept 
    pointOnLine(Px, Py, Ax, Ay, Bx, By) 
    } 

pointOnLine <- function(Px, Py, Ax, Ay, Bx, By) { 

    # This approach based upon comingstorm's answer on 
    # stackoverflow.com/questions/3120357/get-closest-point-to-a-line 
    # Vectorized by Bryan 

    PB <- data.frame(x = Px - Bx, y = Py - By) 
    AB <- data.frame(x = Ax - Bx, y = Ay - By) 
    PB <- as.matrix(PB) 
    AB <- as.matrix(AB) 
    k_raw <- k <- c() 
    for (n in 1:nrow(PB)) { 
     k_raw[n] <- (PB[n,] %*% AB[n,])/(AB[n,] %*% AB[n,]) 
     if (k_raw[n] < 0) { k[n] <- 0 
      } else { if (k_raw[n] > 1) k[n] <- 1 
       else k[n] <- k_raw[n] } 
     } 
    x = (k * Ax + (1 - k)* Bx) 
    y = (k * Ay + (1 - k)* By) 
    ans <- data.frame(x, y) 
    ans 
    } 

# The following proves that pointOnLineNearPoint 
# and pointOnLine work properly and accept vectors 

par(mar = c(4, 4, 4, 4)) # otherwise the plot is slightly distorted 
# and right angles don't appear as right angles 

m <- runif(1, -5, 5) 
b <- runif(1, -20, 20) 
plot(-20:20, -20:20, type = "n", xlab = "x values", ylab = "y values") 
abline(b, m) 

Px <- rnorm(10, 0, 4) 
Py <- rnorm(10, 0, 4) 

res <- pointOnLineNearPoint(Px, Py, m, b) 
points(Px, Py, col = "red") 
segments(Px, Py, res[,1], res[,2], col = "blue") 

##======================================================== 
## 
## Credits: 
## Theory by Paul Bourke http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/ 
## Based in part on C code by Damian Coventry Tuesday, 16 July 2002 
## Based on VBA code by Brandon Crosby 9-6-05 (2 dimensions) 
## With grateful thanks for answering our needs! 
## This is an R (http://www.r-project.org) implementation by Gregoire Thomas 7/11/08 
## 
##======================================================== 

Mã Second Chunk (Plots trình diễn):

set.seed(55) 
np <- 10 # number of data points 
x <- 1:np 
e <- rnorm(np, 0, 60) 
y <- 12 + 5 * x + e 

par(mar = c(4, 4, 4, 4)) # otherwise the plot is slightly distorted 

plot(x, y, main = "Regression minimizes the y-residuals & PCA the normals") 
yx.lm <- lm(y ~ x) 
lines(x, predict(yx.lm), col = "red", lwd = 2) 
segments(x, y, x, fitted(yx.lm), col = "pink") 

# pca "by hand" 
xyNorm <- cbind(x = x - mean(x), y = y - mean(y)) # mean centers 
xyCov <- cov(xyNorm) 
eigenValues <- eigen(xyCov)$values 
eigenVectors <- eigen(xyCov)$vectors 

# Add the first PC by denormalizing back to original coords: 

new.y <- (eigenVectors[2,1]/eigenVectors[1,1] * xyNorm[x]) + mean(y) 
lines(x, new.y, col = "blue", lwd = 2) 

# Now add the normals 

yx2.lm <- lm(new.y ~ x) # zero residuals: already a line 
res <- pointOnLineNearPoint(x, y, yx2.lm$coef[2], yx2.lm$coef[1]) 
points(res[,1], res[,2], col = "blue", pch = 20) # segments should end here 
segments(x, y, res[,1], res[,2], col = "lightblue1") # the normals 
############ CẬP NHẬT

trong lúc 0.123.Tôi tìm thấy gần như chính xác những gì tôi muốn. Nhưng, nó không hoàn toàn hoạt động (rõ ràng là được sử dụng để làm việc). Dưới đây là một mã số đoạn trích từ trang web đó mà âm mưu normals với máy PC đầu tiên phản ánh thông qua một trục dọc:

set.seed(1) 
x <- rnorm(20) 
y <- x + rnorm(20) 
plot(y~x, asp = 1) 
r <- lm(y~x) 
abline(r, col='red') 

r <- princomp(cbind(x,y)) 
b <- r$loadings[2,1]/r$loadings[1,1] 
a <- r$center[2] - b * r$center[1] 
abline(a, b, col = "blue") 
title(main='Appears to use the reflection of PC1') 

u <- r$loadings 
# Projection onto the first axis 
p <- matrix(c(1,0,0,0), nrow=2) 
X <- rbind(x,y) 
X <- r$center + solve(u, p %*% u %*% (X - r$center)) 
segments(x, y, X[1,], X[2,] , col = "lightblue1") 

Và đây là kết quả:

enter image description here

Trả lời

7

Được rồi, tôi sẽ phải trả lời câu hỏi của riêng mình! Sau khi đọc và so sánh thêm về các phương pháp mà mọi người đã đưa lên internet, tôi đã giải quyết được vấn đề. Tôi không chắc chắn tôi có thể nêu rõ những gì tôi "cố định" bởi vì tôi đã trải qua một vài lần lặp lại. Dù sao, đây là cốt truyện và mã (MWE). Các hàm trợ giúp ở cuối cho rõ ràng.

Working Demo

# Comparison of Linear Regression & PCA 
# Generate sample data 

set.seed(39) # gives a decent-looking example 
np <- 10 # number of data points 
x <- -np:np 
e <- rnorm(length(x), 0, 10) 
y <- rnorm(1, 0, 2) * x + 3*rnorm(1, 0, 2) + e 

# Plot the main data & residuals 

plot(x, y, main = "Regression minimizes the y-residuals & PCA the normals", asp = 1) 
yx.lm <- lm(y ~ x) 
lines(x, predict(yx.lm), col = "red", lwd = 2) 
segments(x, y, x, fitted(yx.lm), col = "pink") 

# Now the PCA using built-in functions 
# rotation = loadings = eigenvectors 

r <- prcomp(cbind(x,y), retx = TRUE) 
b <- r$rotation[2,1]/r$rotation[1,1] # gets slope of loading/eigenvector 1 
a <- r$center[2] - b * r$center[1] 
abline(a, b, col = "blue") # Plot 1st PC 

# Plot normals to 1st PC 

X <- pointOnLineNearPoint(x, y, b, a) 
segments(x, y, X[,1], X[,2], col = "lightblue1") 

###### Needed Functions 

pointOnLineNearPoint <- function(Px, Py, slope, intercept) { 
    # Px, Py is the point to test, can be a vector. 
    # slope, intercept is the line to check distance. 

    Ax <- Px-10*diff(range(Px)) 
    Bx <- Px+10*diff(range(Px)) 
    Ay <- Ax * slope + intercept 
    By <- Bx * slope + intercept 
    pointOnLine(Px, Py, Ax, Ay, Bx, By) 
    } 

pointOnLine <- function(Px, Py, Ax, Ay, Bx, By) { 

    # This approach based upon comingstorm's answer on 
    # stackoverflow.com/questions/3120357/get-closest-point-to-a-line 
    # Vectorized by Bryan 

    PB <- data.frame(x = Px - Bx, y = Py - By) 
    AB <- data.frame(x = Ax - Bx, y = Ay - By) 
    PB <- as.matrix(PB) 
    AB <- as.matrix(AB) 
    k_raw <- k <- c() 
    for (n in 1:nrow(PB)) { 
     k_raw[n] <- (PB[n,] %*% AB[n,])/(AB[n,] %*% AB[n,]) 
     if (k_raw[n] < 0) { k[n] <- 0 
      } else { if (k_raw[n] > 1) k[n] <- 1 
       else k[n] <- k_raw[n] } 
     } 
    x = (k * Ax + (1 - k)* Bx) 
    y = (k * Ay + (1 - k)* By) 
    ans <- data.frame(x, y) 
    ans 
    } 
1

Hãy thử thay đổi dòng này của mã của bạn :

res <- pointOnLineNearPoint(x, y, yx2.lm$coef[2], yx2.lm$coef[1]) 

để

res <- pointOnLineNearPoint(x, new.y, yx2.lm$coef[2], yx2.lm$coef[1]) 

Vì vậy, bạn đang gọi đúng các giá trị y.

+0

Ah, tôi có thể không rõ ràng. Các đường màu xanh nhạt phải vuông góc (bình thường) với đường màu xanh, và bắt đầu trên dữ liệu gốc (vòng tròn mở màu đen). Cảm ơn. –

1

Trong Vincent Zoonekynd's code, thay đổi dòng u <- r$loadings để u <- solve(r$loadings). Trong trường hợp thứ hai của solve(), điểm thành phần được dự đoán dọc theo trục chính đầu tiên (tức là ma trận điểm dự đoán với điểm thành phần dự đoán thứ hai được đặt bằng 0) cần được nhân với nghịch đảo của tải/số liệu riêng. Nhân dữ liệu bằng tải cho điểm dự đoán; chia điểm dự đoán cho các lần tải cho dữ liệu. Hy vọng rằng sẽ giúp.

+0

Thú vị. Cảm ơn bạn. Tôi chắc là mã của Vincent đã từng hoạt động. Tôi tự hỏi làm thế nào vấn đề có trong đó. Anh ấy phải đã đăng mã dự thảo nhưng một con số từ mã cuối cùng. –

Các vấn đề liên quan