Here we use a simulated dataset to illustrate the process of Principal Component Analysis (PCA). The original data is three-dimensional, and it was simulated based on a plane described by \(z = 2x+3y\), with some added random noise. As a result, this dataset can be assumed to lie in a two-dimensional space, with two independent vectors defining the surface.

# Load necessary libraries
library(ggplot2)
library(scatterplot3d)

# Set the seed for random number generation
set.seed(123)

# Number of samples (data points) to generate
n_samples <- 20

# Generate n_samples random numbers for x and y in the range [-1, 1] from a uniform distribution
x <- runif(n_samples, -1, 1)
y <- runif(n_samples, -1, 1)

# Define a plane equation z = 2 * x + 3 * y

# z0 represents the ideal values for z without any noise
z0 <- 2 * x + 3 * y 
# Generate random noise with mean 0 and standard deviation 1
z <- z0 + rnorm(n_samples, mean = 0, sd = 1)

# Combine x, y, and z
data_3D <- data.frame(x, y, z)

head(data_3D)
##            x         y          z
## 1 -0.4248450 0.7790786  0.4197223
## 2  0.5766103 0.3856068  2.0920661
## 3 -0.1820462 0.2810136 -0.5470559
## 4  0.7660348 0.9885396  3.7687970
## 5  0.8809346 0.3114116  2.0710647
## 6 -0.9088870 0.4170609 -2.2532845

By visualizing the data, we can observe that the 20 data points are roughly distributed around the plane.

s3d <- scatterplot3d(data_3D$x, data_3D$y, data_3D$z, 
              color = "blue", pch = 16, 
              main = "Original 3D Data", 
              xlab = "X Axis", ylab = "Y Axis", zlab = "Z Axis")

s3d$plane3d(Intercept = 0, x.coef = 2, y.coef = 3, col = "lightblue", lty = 2, lwd = 1.2)

Then, by PCA analysis, the data is reduced to two dimensions.

pca_result <- prcomp(data_3D, center = TRUE, scale. = FALSE)

data_2D <- as.data.frame(pca_result$x[, 1:2])
colnames(data_2D) <- c("PC1", "PC2")

ggplot(data_2D, aes(x = PC1, y = PC2)) +
  geom_point(color = "red", alpha = 2) +
  ggtitle("PCA Reduced 2D Data") +
  xlab("Principal Component 1") +
  ylab("Principal Component 2") +
  theme_minimal()

but…

data0 <- data.frame(x, y, z0)

# Define two vectors (v_x and v_y) that lie on the plane in 3D space.
v_x <- c(-0.1, 0, -0.2)
v_y <- c(0, -0.1, -0.3)

# Create a transformation matrix by the two vectors v_x and v_y
transformation_matrix <- cbind(v_x, v_y)[1:2, ]

# Apply the linear transformation to the original data
transformed_data <- data.frame(
  x_transformed = data0$x * v_x[1] + data0$y * v_y[1],
  y_transformed = data0$x * v_x[2] + data0$y * v_y[2] )

ggplot(transformed_data, aes(x = x_transformed, y = y_transformed)) +
  geom_point(color = "blue") +
  labs(title = "2D Visualization of Data after Linear Transformation", x = "Transformed X", y = "Transformed Y") +
  theme_minimal()

they are quite differnet!

We perform the PCA step by step.

0.1 1. standardization

\[ Z_{ij} = \frac{X_{ij} - \mu_{j}}{\sigma_{j}} \]

where \(i\) is the order of data, \(j\) is the dimension.

data_standardized <- scale(data_3D)
head(data_standardized)
##               x         y           z
## [1,] -0.8397295 1.0912011 -0.07826596
## [2,]  0.7576354 0.4049295  0.60138989
## [3,] -0.4524548 0.2225039 -0.47117346
## [4,]  1.0597758 1.4565312  1.28282873
## [5,]  1.2430460 0.2755223  0.59285473
## [6,] -1.6117978 0.4597900 -1.16460041

0.2 2. covariance matrix

\[ Cov = \frac{1}{n-1} X^T X \]

cov_matrix <- cov(data_standardized)

cov_matrix
##           x         y         z
## x 1.0000000 0.2866221 0.7536307
## y 0.2866221 1.0000000 0.7658283
## z 0.7536307 0.7658283 1.0000000

0.3 3. eigenvalue and eigenvector of the cov matrix

\[ Cv = \lambda v \]

eigen_result <- eigen(cov_matrix)
eigenvalues <- eigen_result$values
eigenvectors <- eigen_result$vectors

eigen_result
## eigen() decomposition
## $values
## [1] 2.22726773 0.71342093 0.05931134
## 
## $vectors
##            [,1]         [,2]       [,3]
## [1,] -0.5293650  0.713684741  0.4587230
## [2,] -0.5346723 -0.700449256  0.4727540
## [3,] -0.6587095 -0.004992959 -0.7523808

0.4 4. choose dimensions

sorted_indices <- order(eigenvalues, decreasing = TRUE)
top_2_eigenvectors <- eigenvectors[, sorted_indices[1:2]]

sorted_indices
## [1] 1 2 3

Here are the two eigenvectors we selected (red and green), with a set of independent vectors that lie on the plane (dark red \((1,0,2)\) and dark green \((0.6,-0.5,-0.3)\)).

s3d <- scatterplot3d(data_3D$x, data_3D$y, data_3D$z,
                     color = "blue", pch = 16, 
                     main = "3D Data with PCA Components",
                     xlab = "X Axis", ylab = "Y Axis", zlab = "Z Axis")

s3d$plane3d(Intercept = 0, x.coef = 2, y.coef = 3, col = "lightblue", lty = 2, lwd = 1.2)

s3d$points3d(c(0, top_2_eigenvectors[1, 1]), 
             c(0, top_2_eigenvectors[2, 1]), 
             c(0, top_2_eigenvectors[3, 1]), 
             col = "red", type = "l", lwd = 2)

s3d$points3d(c(0, top_2_eigenvectors[1, 2]), 
             c(0, top_2_eigenvectors[2, 2]), 
             c(0, top_2_eigenvectors[3, 2]), 
             col = "green", type = "l", lwd = 2)

s3d$points3d(c(0, 0.6), 
             c(0, 0), 
             c(0, 1.2), 
             col = "darkred", type = "l", lwd = 5)

s3d$points3d(c(0, 0.6), 
             c(0, -0.5), 
             c(0, 0.3), 
             col = "darkgreen", type = "l", lwd = 5)

s3d$points3d(data_3D$x, data_3D$y, data_3D$z, col = "blue", pch = 16)

0.5 5. Mapping the data to the new space

pca_result <- data_standardized %*% top_2_eigenvectors

head(pca_result)
##             [,1]       [,2]
## [1,] -0.08735706 -1.3632424
## [2,] -1.01371148  0.2540776
## [3,]  0.43091352 -0.4764102
## [4,] -2.18478661 -0.2702855
## [5,] -1.19585824  0.6911934
## [6,]  1.37452571 -1.4665602
data_2D <- as.data.frame(pca_result[, 1:2])
colnames(data_2D) <- c("PC1", "PC2")

ggplot(data_2D, aes(x = PC1, y = PC2)) +
  geom_point(color = "red", alpha = 5) +
  ggtitle("PCA Reduced 2D Data") +
  xlab("Principal Component 1") +
  ylab("Principal Component 2") +
  theme_minimal()

Here we could explore more about the sample size and data noise (which is the variance \(\sigma\) here).

pca_test<-function(n,sigma){
  n_samples <- n
  x <- runif(n_samples, -1, 1)  
  y <- runif(n_samples, -1, 1)
  z0 <- 2 * x + 3 * y 
  z <- z0 + rnorm(n_samples, mean = 0, sd = sigma)
  data_3D <- data.frame(x, y, z)
  
  v_x <- c(-0.1, 0, -0.2)
  v_y <- c(0, -0.1, -0.3)
  
  data_standardized <- scale(data_3D)
  cov_matrix <- cov(data_standardized)
  eigen_result <- eigen(cov_matrix)
  eigenvalues <- eigen_result$values
  eigenvectors <- eigen_result$vectors
  
  sorted_indices <- order(eigenvalues, decreasing = TRUE)
  top_2_eigenvectors <- eigenvectors[, sorted_indices[1:2]]
  
  s3d <- scatterplot3d(data_3D$x, data_3D$y, data_3D$z,
                     color = "blue", pch = 16, 
                     main = "3D Data with PCA Components",
                     xlab = "X Axis", ylab = "Y Axis", zlab = "Z Axis", type = "n")

s3d$plane3d(Intercept = 0, x.coef = 2, y.coef = 3, col = "lightblue", lty = 2, lwd = 1.2)

s3d$points3d(c(0, top_2_eigenvectors[1, 1]), 
             c(0, top_2_eigenvectors[2, 1]), 
             c(0, top_2_eigenvectors[3, 1]), 
             col = "red", type = "l", lwd = 2)

s3d$points3d(c(0, top_2_eigenvectors[1, 2]), 
             c(0, top_2_eigenvectors[2, 2]), 
             c(0, top_2_eigenvectors[3, 2]), 
             col = "green", type = "l", lwd = 2)

s3d$points3d(c(0, 0.6), 
             c(0, 0), 
             c(0, 1.2), 
             col = "darkred", type = "l", lwd = 5)

s3d$points3d(c(0, 0.6), 
             c(0, -0.5), 
             c(0, -0.3), 
             col = "darkgreen", type = "l", lwd = 5)

s3d$points3d(data_3D$x, data_3D$y, data_3D$z, col = "blue", pch = 16, type = "n")
}
#set.seed(2)
pca_test(10000,1)

#set.seed(2)
pca_test(10000,1)

pca_test(10000,1)

pca_test(10,1)

pca_test(10,1)

pca_test(10,1)