2011-12-10 17 views
9

मैं Cerebral Mastication ब्लॉग से प्रेरित प्रतिगमन और पीसीए की तुलना करने के लिए एक विधि को सही करने की कोशिश कर रहा हूं जिस पर SO पर एक अलग कोण से भी चर्चा की गई है। इससे पहले कि मैं भूल जाऊं, जेडी लांग और जोश उलरिच के बहुत सारे कारणों के लिए बहुत धन्यवाद। मैं अगले सत्र में एक कोर्स में इसका उपयोग करने जा रहा हूं। क्षमा करें यह लंबा है!रिग्रेशन और पीसीए की दृश्य तुलना

अद्यतन: मुझे एक अलग दृष्टिकोण मिला जो लगभग काम करता है (अगर आप कर सकते हैं तो इसे ठीक करें!)। मैंने इसे नीचे पोस्ट किया। मैं उससे ज्यादा स्मार्ट और छोटे दृष्टिकोण के साथ आने में सक्षम था!

मैं मूल रूप से पिछली योजनाओं को एक बिंदु तक चलाता हूं: यादृच्छिक डेटा उत्पन्न करें, सर्वोत्तम फिट की रेखा को समझें, अवशिष्ट बनाएं। यह नीचे दूसरे कोड खंड में दिखाया गया है। लेकिन मैंने चारों ओर खोद दिया और एक यादृच्छिक बिंदु (इस मामले में डेटा पॉइंट) के माध्यम से लाइनों को सामान्य रेखाओं को आकर्षित करने के लिए कुछ फ़ंक्शन लिखे। मुझे लगता है कि ये काम ठीक है, और वे सबूत के साथ पहले कोड चंक में दिखाए जाते हैं।

अब, दूसरा कोड चंक @JDLong के समान प्रवाह का उपयोग करके कार्रवाई में पूरी चीज दिखाता है और मैं परिणामी साजिश की एक छवि जोड़ रहा हूं। काला, लाल रंग का डेटा गुलाबी अवशेषों के साथ प्रतिगमन है, नीला पहला पीसी है और हल्का नीला मानक होना चाहिए, लेकिन जाहिर है कि वे नहीं हैं। इन नियमों को आकर्षित करने वाले फर्स्ट कोड चंक में कार्य ठीक लगते हैं, लेकिन प्रदर्शन के साथ कुछ सही नहीं है: मुझे लगता है कि मुझे कुछ गलत समझना होगा या गलत मूल्यों को पार करना होगा। मेरे आदर्श क्षैतिज में आते हैं, जो एक उपयोगी सुराग की तरह लगता है (लेकिन अब तक, मेरे लिए नहीं)। क्या कोई देख सकता है कि यहां क्या गलत है?

धन्यवाद, यह थोड़ी देर के लिए मुझे परेशान किया गया है ... Plot showing problem

पहले कोड हिस्सा (कार्य नॉर्मल्स ड्रा करने के लिए और सबूत वे काम):

##### 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 
## 
##======================================================== 

दूसरा कोड हिस्सा (भूखंड प्रदर्शन):

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 
############ अद्यतन

पर अधिकमुझे लगभग वही मिला जो मैं चाहता था। लेकिन, यह काफी काम नहीं करता है (स्पष्ट रूप से काम करने के लिए प्रयोग किया जाता है)।

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") 

और यहाँ परिणाम है::

enter image description here

उत्तर

7

ठीक है, मुझे अपने प्रश्न का उत्तर देना होगा! लोगों ने इंटरनेट पर रखे तरीकों की और पढ़ने और तुलना करने के बाद, मैंने समस्या हल कर ली है। मुझे यकीन नहीं है कि मैं स्पष्ट रूप से बता सकता हूं कि मैंने "तय" किया क्योंकि मैं काफी कुछ पुनरावृत्तियों से गुजर चुका हूं। वैसे भी, यहां साजिश और कोड (MWE) है। सहायकता के लिए सहायक कार्य अंत में हैं।

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

अपने कोड की इस पंक्ति बदलने का प्रयास करें यहाँ कि साइट है जो भूखंडों पहले पीसी के लिए normals एक ऊर्ध्वाधर अक्ष के माध्यम से परिलक्षित से एक कोड अंश है :

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

को

तो आप सही वाई मानों को बुला रहे हैं।

+0

आह, मैं स्पष्ट नहीं हो सकता था। हल्की नीली रेखाएं नीली रेखा के लिए लंबवत (सामान्य) होनी चाहिए, और मूल डेटा (काला खुली मंडल) से शुरू होनी चाहिए। धन्यवाद। –

1

Vincent Zoonekynd's code में, u <- solve(r$loadings) करने के लिए लाइन u <- r$loadings बदल जाते हैं। solve() के दूसरे उदाहरण में, पहले प्रिंसिपल अक्ष के साथ अनुमानित घटक स्कोर (यानी, अनुमानित स्कोर के मैट्रिक्स को दूसरे अनुमानित घटक स्कोर शून्य पर सेट किया गया है) लोडिंग/ईजिनवेक्टर के उलटा द्वारा गुणा करने की आवश्यकता है। लोडिंग द्वारा डेटा गुणा करने से अनुमानित स्कोर मिलते हैं; लोडिंग द्वारा अनुमानित स्कोर को विभाजित करने से डेटा मिलता है। उम्मीद है की वो मदद करदे।

+0

दिलचस्प। धन्यवाद। मुझे यकीन है कि विन्सेंट का कोड काम करने के लिए प्रयोग किया जाता है। मुझे आश्चर्य है कि समस्या में कैसे समस्या आई। उन्होंने मसौदा कोड पोस्ट किया होगा लेकिन अंतिम कोड से एक आंकड़ा होना चाहिए। –

संबंधित मुद्दे