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 ...
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ả:
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. –